% FEGEOMETRY Creates an fegeometry object that can be used for finite
% element analysis
%
% gm = fegeometry(INPUT) creates an fegeometry object based on the
% specified INPUT. The INPUT can be:
%       - STL or STEP file
%       - function handle
%       - geometry description matrix
%       - triangulation object
%       - delaunayTriangulation object
%       - DiscreteGeometry object
%       - AnalyticGeometry object
%
% gm = fegeometry(Nodes,Elements) creates an fegeometry object from
% a mesh defined by Nodes and Elements.
%
% Here, Nodes is an Nn-by-2 or Nn-by-3 matrix, and Elements is an
% Ne-by-Nnpe matrix. Nn is the number of nodes in the mesh. Ne is the
% number of elements in the mesh. Nnpe is the number of nodes per element,
% which is 3 or 6 for a planar or surface mesh, and 4 or 10 for a volume
% mesh. If the input is a planar or volume mesh, it is assigned to the Mesh
% property.
%
% gm = fegeometry(Nodes,Elements,RegionId) constructs a multidomain
% geometric model from a mesh defined by Nodes and Elements. Domain
% information is specified by a vector RegionId. The length of this vector
% equals the number of elements in the mesh. RegionId(i) is the ID of the
% geometric region associated with element i.
%
% fegeometry 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
%    Vertices - The coordinates of the vertices in the geometry
%    Mesh - FEMesh object containing a description of the finite element mesh
%
%
% fegeometry methods
%
%    generateMesh - generates a mesh for fegeometry
%    pdemesh - plots the mesh of the geometry
%    pdegplot - plots fegeometry
%    triangulation - return a triangulation object from the fegeometry
%
%  Geometry modification:
%    addVertex
%    addFace
%    addCell
%    addVoid
%    extrude
%    mergeCells
%    rotate
%    scale
%    translate
%
%  Geometry queries:
%    cellFaces
%    facesAttachedToEdges
%    nearestFace
%    cellEdges
%    faceEdges
%    nearestEdge
%
%
% Example 1: Create a geometry from an STL file.
% gm = fegeometry('ForearmLink.stl');
%
% Example 2: Create a geometry from a function handle.
% gm = fegeometry(@cardg);
%
% Example 3: Create a geometry from a geometry description matrix.
% g=[3 4 0 1 1 0 0 0 1.0 1.0];
% sf = 'S1';
% ns = 'S1';
% gm = fegeometry(decsg(g', sf, ns'));
%
% Example 4: Create a geometry from a DiscreteGeometry object.
% gm = fegeometry(multicuboid(1,1,1));
%
% Example 5: Create a geometry from an AnalyticGeometry object, generated
% by a call to geometryFromEdges.
% agm = geometryFromEdges(@lshapeg);
% gm = fegeometry(agm);
%
% Example 6: Create a geometry from node and element descriptions.
% load tetmesh;
% gm = fegeometry(X,tet);

% Copyright 2023 The MathWorks, Inc.

classdef fegeometry

    properties (Dependent,SetAccess=private)
        NumFaces
        NumEdges
        NumVertices
        NumCells
        Vertices
    end

    properties (Dependent)
        Mesh
    end

    properties (SetAccess=private,Hidden)
        GeomRep
        MeshRep
    end

    properties (Dependent,SetAccess=private,Hidden)
        Dim
    end

    methods
        function self = fegeometry(input,elements,regions,options)
            arguments
                input {fegeometry.mustBeValidInput} = []
                elements (:,:) double {mustBeNonsparse,mustBeInteger,fegeometry.mustMatchElementSize} = []
                regions (1,:) double {mustBeNonsparse,mustBeInteger} = []
                options.MaxRelativeDeviation (1,1) double {mustBeInRange(options.MaxRelativeDeviation,0.1,10)}
                options.AllowSelfIntersections (1,1) logical {mustBeNumericOrLogical} = true
                options.FeatureAngle (1,1) double {mustBeInRange(options.FeatureAngle,10,90)} = 44
            end

            try
                if nargin == 0 % empty constructor for object arrays
                    self.GeomRep = [];
                elseif nargin > 1
                    nodes = input';
                    elems = elements';
                    checkMeshValidity(nodes,elems,regions);
                    if nargin == 2
                        [self.GeomRep, self.Mesh] = geomFromMeshInternal(nodes,elems);
                    elseif isempty(regions)
                        error(message('pde:pdeModel:InvalidSubdomainId'));
                    else
                        [self.GeomRep, self.Mesh] = geomFromMeshInternal(nodes,elems,regions);
                    end
                elseif isempty(input)
                    self = fegeometry.empty;
                elseif isa(input,'char') || (isscalar(input) && isstring(input))
                    importopts = namedargs2cell(options);
                    self.GeomRep = importGeometry(input,importopts{:});
                elseif isa(input,'function_handle') || isnumeric(input)
                    self.GeomRep = self.copy(geometryFromEdges(input));
                elseif isa(input,'triangulation') % includes DelaunayTriangulation
                    nodes = input.Points';
                    elems = input.ConnectivityList';
                    checkMeshValidity(nodes,elems,regions);
                    [self.GeomRep, self.Mesh] = geomFromMeshInternal(nodes,elems,regions,[],options.FeatureAngle);
                    if ~options.AllowSelfIntersections
                        checkForSelfIntersections(self.GeomRep);
                    end
                else % pde.DiscreteGeometry, pde.AnalyticGeometry
                    self.GeomRep = self.copy(input);
                    if ~options.AllowSelfIntersections && isa(input,"pde.DiscreteGeometry")
                        checkForSelfIntersections(self.GeomRep);
                    end
                end
            catch ex
                throwAsCaller(ex);
            end
        end

        % Dependent properties
        function numfaces = get.NumFaces(self)
            numfaces = self.GeomRep.NumFaces;
        end
        function numedges = get.NumEdges(self)
            numedges = self.GeomRep.NumEdges;
        end
        function numvertices = get.NumVertices(self)
            numvertices = self.GeomRep.NumVertices;
        end
        function numcells = get.NumCells(self)
            numcells = self.GeomRep.NumCells;
        end
        function vertices = get.Vertices(self)
            vertices = self.GeomRep.Vertices;
        end

        function dim = get.Dim(self)
            if isempty(self.GeomRep)
                dim = [];
            elseif self.NumCells < 1
                dim = 2;
            else
                dim = 3;
            end
        end

        function mesh = get.Mesh(self)
            if isempty(self.MeshRep)
                mesh = [];
            else
                meshProps = ["Nodes","Elements","MaxElementSize","MinElementSize",...
                    "MeshGradation","GeometricOrder","MeshAssociation"];
                z = cellfun(@(k) self.MeshRep.(k),meshProps,'UniformOutput',false);
                mesh = pde.FEMesh(z{:});
            end
        end

        function self = set.Mesh(self,mesh)
            if isempty(mesh)
                self.MeshRep = [];
            else
                self.MeshRep = struct();
                meshProps = ["Nodes","Elements","MaxElementSize","MinElementSize",...
                    "MeshGradation","GeometricOrder","MeshAssociation"];
                for k = meshProps
                    self.MeshRep.(k) = mesh.(k);
                end
            end
        end

        % Plot and mesh methods
        function self = generateMesh(self,varargin)
            if ~isempty(self)
                tmp = createpde();
                tmp.Geometry = self.GeomRep;
                self.Mesh = generateMesh(tmp,varargin{:});
            end
            if nargout < 1
                warning(message('pde:fegeometry:NoOutput'));
            end
        end

        function h = pdegplot(self,varargin)
            if ~isempty(self)
                hh = pdegplot(self.GeomRep,varargin{:});
                if nargout > 0
                    h = hh;
                end
            end
        end

        function h = pdemesh(self,varargin)
            if ~isempty(self)
                hh = pdemesh(self.Mesh,varargin{:});
                if nargout > 0
                    h = hh;
                end
            end
        end

        % Geometry modification
        function [self,vid] = addVertex(self,options)
            arguments
                self (1,1) fegeometry
                options.Coordinates (:,:) {mustBeNumeric,mustBeFinite}
            end
            if isempty(self.GeomRep)
                error(message("pde:fegeometry:ModifyEmpty"));
            end
            if isfield(options,"Coordinates")
                options.Coordinates = double(options.Coordinates);
            end
            if isa(self.GeomRep,'pde.AnalyticGeometry')
                gm = convertAnalyticToDiscrete(self.GeomRep);
            else
                gm = self.copy(self.GeomRep);
            end
            args = namedargs2cell(options);
            try
                vid = addVertex(gm,args{:});
            catch ME
                throwAsCaller(ME);
            end
            self.Mesh = [];
            self.GeomRep = gm;
            if nargout < 1
                warning(message('pde:fegeometry:NoOutput'));
            end
        end

        function [self,fid] = addFace(self,edges)
            arguments
                self (1,1) fegeometry
                edges (1,:) % integer vector or cell array of integer vectors
            end
            if isempty(self.GeomRep)
                error(message("pde:fegeometry:ModifyEmpty"));
            end
            gm = self.copy(self.GeomRep);
            try
                [self.GeomRep,fid] = addFace(gm,edges);
            catch ME
                throwAsCaller(ME);
            end
            self.Mesh = [];
            if nargout < 1
                warning(message('pde:fegeometry:NoOutput'));
            end
        end

        function self = addCell(self,geom2)
            arguments
                self (1,1) fegeometry
                geom2 (1,1) fegeometry
            end
            if ~isa(self.GeomRep,'pde.DiscreteGeometry') || ...
                    ~isa(geom2.GeomRep,'pde.DiscreteGeometry')
                error(message('pde:pdeGeometryClass:GeometryMustBe3D'));
            end
            gm = self.copy(self.GeomRep);
            gm2 = self.copy(geom2.GeomRep);
            try
                self.GeomRep = addCell(gm,gm2);
            catch ME
                throwAsCaller(ME);
            end
            self.Mesh = [];
            if nargout < 1
                warning(message('pde:fegeometry:NoOutput'));
            end
        end

        function self = addVoid(self,geom2)
            arguments
                self (1,1) fegeometry
                geom2 (1,1) fegeometry
            end
            if ~isa(self.GeomRep,'pde.DiscreteGeometry') || ...
                    ~isa(geom2.GeomRep,'pde.DiscreteGeometry')
                error(message('pde:pdeGeometryClass:GeometryMustBe3D'));
            end
            gm = self.copy(self.GeomRep);
            gm2 = self.copy(geom2.GeomRep);
            try
                self.GeomRep = addVoid(gm,gm2);
            catch ME
                throwAsCaller(ME);
            end
            self.Mesh = [];
            if nargout < 1
                warning(message('pde:fegeometry:NoOutput'));
            end
        end

        function self = mergeCells(self,cells)
            arguments
                self (1,1) fegeometry
                cells (1,:) {mustBeNumeric,mustBePositive,mustBeInteger}
            end
            if isempty(self.GeomRep)
                error(message("pde:fegeometry:ModifyEmpty"));
            end
            if self.Dim ~= 3
                error(message('pde:pdeGeometryClass:GeometryMustBe3D'));
            end
            if length(cells) ~= 2
                error(message('pde:pdeGeometryClass:CellsMustBeVectorOfTwo'));
            end
            if cells(1) == cells(2)
                error(message('pde:pdeGeometryClass:CellsMustBeUnique'));
            end
            gm = self.copy(self.GeomRep);
            try
                self.GeomRep = mergeCells(gm,cells);
            catch ME
                throwAsCaller(ME);
            end
            self.Mesh = [];
            if nargout < 1
                warning(message('pde:fegeometry:NoOutput'));
            end
        end

        function self = extrude(self,HeightOrFaceID,height)
            arguments
                self (1,1) fegeometry
                HeightOrFaceID (1,:) {mustBeNumeric,mustBePositive,mustBeFinite}
                height (1,:) {mustBeNumeric,mustBePositive,mustBeFinite} = []
            end
            if isempty(self.GeomRep)
                error(message("pde:fegeometry:ModifyEmpty"));
            end
            gm = self.copy(self.GeomRep);
            HeightOrFaceID = double(HeightOrFaceID);
            height = double(height);
            try
                if nargin > 2
                    self.GeomRep = extrude(gm,HeightOrFaceID,height);
                else
                    self.GeomRep = extrude(gm,HeightOrFaceID);
                end
            catch ME
                throwAsCaller(ME);
            end
            self.Mesh = [];
            if nargout < 1
                warning(message('pde:fegeometry:NoOutput'));
            end
        end

        function self = rotate(self,theta,refpoint1,refpoint2)
            arguments
                self (1,1) fegeometry
                theta (1,1) {mustBeNumeric}
                refpoint1 (1,:) {mustBeNumeric} = []
                refpoint2( 1,:) {mustBeNumeric} = []
            end
            if ~isempty(self.GeomRep)
                gm = self.copy(self.GeomRep);
                theta = double(theta);
                refpoint1 = double(refpoint1);
                refpoint2 = double(refpoint2);
                try
                    if ~isempty(refpoint2)
                        gm = rotate(gm,theta,refpoint1,refpoint2);
                    elseif ~isempty(refpoint1)
                        gm = rotate(gm,theta,refpoint1);
                    else
                        gm = rotate(gm,theta);
                    end
                catch ME
                    throwAsCaller(ME);
                end
                self.Mesh = [];
                self.GeomRep = gm;
                if nargout < 1
                    warning(message('pde:fegeometry:NoOutput'));
                end
            end
        end

        function self = scale(self,s,refpoint)
            arguments
                self (1,1) fegeometry
                s (1,:) {mustBeNumeric}
                refpoint (1,:) {mustBeNumeric} = []
            end
            if ~isempty(self.GeomRep)
                gm = self.copy(self.GeomRep);
                s = double(s);
                refpoint = double(refpoint);
                try
                    if isempty(refpoint)
                        gm = scale(gm,s);
                    else
                        gm = scale(gm,s,refpoint);
                    end
                catch ME
                    throwAsCaller(ME);
                end
                self.Mesh = [];
                self.GeomRep = gm;
                if nargout < 1
                    warning(message('pde:fegeometry:NoOutput'));
                end
            end
        end

        function self = translate(self,s)
            arguments
                self (1,1) fegeometry
                s (1,:) {mustBeNumeric}
            end
            if ~isempty(self.GeomRep)
                gm = self.copy(self.GeomRep);
                s = double(s);
                try
                    gm = translate(gm,s);
                catch ME
                    throwAsCaller(ME);
                end
                self.Mesh = [];
                self.GeomRep = gm;
                if nargout < 1
                    warning(message('pde:fegeometry:NoOutput'));
                end
            end
        end

        % Queries
        function v = cellFaces(self,RegionID,FilterType)
            arguments
                self (1,1) fegeometry
                RegionID (1,:) double {mustBePositive,mustBeInteger}
                FilterType (1,1) string {mustBeMember(FilterType,...
                    ["all","internal","external"])} = "all"
            end
            try
                v = cellFaces(self.GeomRep,RegionID,FilterType);
            catch ME
                throwAsCaller(ME);
            end
        end

        function v = facesAttachedToEdges(self,RegionID,FilterType)
            arguments
                self (1,1) fegeometry
                RegionID (1,:) double {mustBePositive,mustBeInteger}
                FilterType (1,1) string {mustBeMember(FilterType,...
                    ["all","internal","external"])} = "all"
            end
            try
                v = facesAttachedToEdges(self.GeomRep,RegionID,FilterType);
            catch ME
                throwAsCaller(ME);
            end
        end

        function v = nearestFace(self,coords)
            arguments
                self (1,1) fegeometry
                coords (:,:) double {mustBeFinite}
            end
            try
                v = nearestFace(self.GeomRep,coords);
            catch ME
                throwAsCaller(ME);
            end
        end

        function v = cellEdges(self,RegionID,FilterType)
            arguments
                self (1,1) fegeometry
                RegionID (1,:) double {mustBePositive,mustBeInteger}
                FilterType (1,1) string {mustBeMember(FilterType,...
                    ["all","internal","external"])} = "all"
            end
            try
                v = cellEdges(self.GeomRep,RegionID,FilterType);
            catch ME
                throwAsCaller(ME);
            end
        end

        function v = faceEdges(self,RegionID,FilterType)
            arguments
                self (1,1) fegeometry
                RegionID (1,:) double {mustBePositive,mustBeInteger}
                FilterType (1,1) string {mustBeMember(FilterType,...
                    ["all","internal","external"])} = "all"
            end
            try
                v = faceEdges(self.GeomRep,RegionID,FilterType);
            catch ME
                throwAsCaller(ME);
            end
        end

        function v = nearestEdge(self,coords)
            arguments
                self (1,1) fegeometry
                coords (:,:) double {mustBeFinite}
            end
            try
                v = nearestEdge(self.GeomRep,coords);
            catch ME
                throwAsCaller(ME);
            end
        end

        function tri = triangulation(self)
            % TRIANGULATION - return a Triangulation object from fegeometry
            arguments
                self (1,1) fegeometry
            end
            % get discrete geometry
            if isempty(self.GeomRep)
                error(message('pde:fegeometry:GeoEmpty'));
            end
            if isa(self.GeomRep,'pde.AnalyticGeometry')
                gm = convertAnalyticToDiscrete(self.GeomRep);
            else
                gm = self.GeomRep;
            end
            ws = warning('off','MATLAB:structOnObject');
            restoreWarning = onCleanup(@() warning(ws));
            gstruct = struct(gm);
            % check if multidomain
            if (gstruct.IsTwoD && gstruct.NumFaces > 1) || ...
                    (~gstruct.IsTwoD && gstruct.NumCells > 1)
                error(message('pde:fegeometry:TriMultidomain'));
            end
            % convert triangulation
            gtri = gstruct.Triangulation;
            gassoc = gstruct.Associativity;
            if gstruct.IsTwoD
                % remove zero points from triangulation
                gtri.Points = gtri.Points(:,1:2);
                tri = triangulation(gtri.Triangles,gtri.Points);
            else
                tri = cellTriangulation(gtri,gassoc,1);
            end        
        end
    end

    methods(Static,Access=private)
        function newObj = copy(obj)
            objByteArray = getByteStreamFromArray(obj);
            newObj = getArrayFromByteStream(objByteArray);
        end

        function mustMatchElementSize(x)
            if ~isempty(x) && ~ismember(size(x,2),[3,4,6,10])
                error(message('pde:fegeometry:WrongElementSize'));
            end
        end

        function mustBeValidInput(x)
            if ~(ischar(x) && isvector(x)) && ...
                    ~(isscalar(x) && isstring(x)) && ...
                    ~isnumeric(x) && ~isa(x,'function_handle') && ...
                    ~isempty(x) && ...
                    ~(isscalar(x) && isa(x,'pde.DiscreteGeometry')) && ...
                    ~(isscalar(x) && isa(x,'pde.AnalyticGeometry')) && ...
                    ~isa(x,'triangulation')
                error(message('pde:fegeometry:InvalidInput'));
            end
        end
    end
end

function tri = cellTriangulation(tr,assoc,domains)
% Extract a single cell from a multidomain geometry.

% collect faces and their sense codes for each entry in domains
faceIndices = indexCounts(assoc.CellFaceCounts,domains);
faces = assoc.CellFaces(faceIndices,1);
facecodes = assoc.CellFaces(faceIndices,2);

% classify faces
bwdFaces = faces(facecodes == 0) + 1;
fwdFaces = faces(facecodes > 0) + 1;

% get facet lists for both sets of faces
fwdFacets = indexCounts(assoc.FaceFacetCounts,fwdFaces);
bwdFacets = indexCounts(assoc.FaceFacetCounts,bwdFaces);
fwdFacets = assoc.FaceFacets(fwdFacets) + 1;
bwdFacets = assoc.FaceFacets(bwdFacets) + 1;

% flip backwards facets
cl1 = tr.Triangles(fwdFacets,:);
cl2 = tr.Triangles(bwdFacets,:);
cl = [cl1; cl2(:,[2 1 3])];
pts = tr.Points;

% remove any unreferenced vertices
ind = unique(cl(:));
pp(ind) = 1:length(ind);
cl = pp(cl);
pts = pts(ind,:);
tri = triangulation(cl,pts);
end

function out = indexCounts(counts,vals)
ptr = cumsum([1;counts]);
ranges = arrayfun(@(x)ptr(x):ptr(x+1)-1,vals,'UniformOutput',false);
out = [ranges{:}];
end

function gdisc = convertAnalyticToDiscrete(grep)
mdl = createpde;
facets = facetAnalyticGeometry(mdl,grep,0);
gdisc = analyticToDiscrete(facets);
end
