function msh = generateMesh(self, varargin)    
% generateMesh            - Generate a mesh from the geometry        
%    MSH = generateMesh(PDEM) Generates a mesh object from the geometry stored        
%    in the Geometry property. The mesh created is assigned to the Mesh 
%    property of this class and a handle to the mesh is returned in MSH.
%
%    MSH = generateMesh(...,'Name1',Value1, 'Name2',Value2, ...) 
%    uses name/value pairs to override the default meshing values for 
%    'Name1', 'Name2',... Supported options for generateMesh include:
%
%    Property        Value/{Default}         Description
%    -----------------------------------------------------------
%    Hmax            numeric {estimate}    Target maximum edge size
%    Hmin            numeric {estimate}    Target minimum edge size
%    Hgrad           numeric {1.5}         Element gradation control
%    Hface           numeric {none}        Target size on selected faces
%    Hedge           numeric {none}        Target size on selected edges
%    Hvertex         numeric {none}        Target size at selected vertices
%    GeometricOrder  linear|{quadratic}    Element geometric order 
%                                          
%
%    Note: If Hmax is inadvertently set too low the mesh generation may 
%    take a long time to complete. Use Ctrl-C to interrupt the meshing.
%    The mesh generator may take additional time to respond to the
%    interrupt.
%
%    Example 1:  Import a 3D geometry, generate the default mesh and plot.
%       numberOfPDEs = 1;
%       pdem = createpde(numberOfPDEs);
%       gm = importGeometry(pdem,'Block.stl')
%       msh = generateMesh(pdem)
%       pdemesh(msh);
%
%    Example 2:  Import a 3D geometry, generate the locally refined mesh and plot.
%       pdem = createpde();
%       gm = importGeometry(pdem,'Blade.stl')
%       msh = generateMesh(pdem,'Hedge',{[20 21],0.002})
%       pdeplot3D(msh);
%
% See also pdeplot, pdeplot3D, pdemesh, pde.FEMesh

% Copyright 2014-2021 The MathWorks, Inc.

    if nargin > 1
        [varargin{:}] = convertStringsToChars(varargin{:});
    end

    geom = self.Geometry;
    if isempty(geom)
       msh = [];
       return;
    end    
    
    
    p = inputParser;
    addParameter(p,'Hmax',0,@pde.PDEModel.isValidHmax);
    addParameter(p,'Hmin',0,@pde.PDEModel.isValidHmin);       
    addParameter(p,'Hgrad',0,@pde.PDEModel.isValidHgrad);
    addParameter(p,'GeometricOrder','quadratic', @pde.PDEModel.isValidGeomOrder);        
    addParameter(p,'Hface',{},@pde.PDEModel.isValidHface);
    addParameter(p,'Hedge',{},@pde.PDEModel.isValidHedge);
    addParameter(p,'Hvertex',{},@pde.PDEModel.isValidHvertex);

    if self.IsTwoD && isa(geom, 'pde.AnalyticGeometry')
        addParameter(p,'MesherVersion','');   % No longer used
        addParameter(p,'Jiggle','');          % No longer used
        addParameter(p,'JiggleIter','');      % No longer used      
    end
    parse(p,varargin{:});
    Hmax = p.Results.Hmax;
    Hmin = p.Results.Hmin;
    Hgrad = p.Results.Hgrad;
    Hface = p.Results.Hface;
    Hedge = p.Results.Hedge;
    Hvertex = p.Results.Hvertex;
    geomOrder = convertStringsToChars(p.Results.GeometricOrder);

    % sanity checks
    if Hmax ~= 0 && Hmin ~= 0 && Hmin > Hmax         
        error(message('pde:pdeModel:invalidHmaxHmin'));  
    end
    nc = numel(geomOrder);        
    if (strncmpi(geomOrder,'linear',nc))
        geomOrder = 'linear';
    else
        geomOrder = 'quadratic';
    end

    % strip and check local refinement data
    Hface = getValidLocalSizes(Hface,geom.NumFaces,'pde:pdeGeometryClass:InvalidFaceIndex','pde:pdeModel:invalidHface');
    Hedge = getValidLocalSizes(Hedge,geom.NumEdges,'pde:pdeGeometryClass:InvalidEdgeIndex','pde:pdeModel:invalidHedge');
    Hvertex = getValidLocalSizes(Hvertex,geom.NumVertices,'pde:pdeGeometryClass:InvalidVertexIndex','pde:pdeModel:invalidHvertex');
    
    if self.IsTwoD && isa(geom, 'pde.AnalyticGeometry')                                         
        try
            msh = generateMeshFacettedAnalytic(self,Hmax,Hmin,Hgrad,Hface,Hedge,Hvertex,geomOrder);
        catch ex
            throwAsCaller(ex);
         end   
        self.Mesh = msh;
        addlistener(msh,'ObjectBeingDestroyed',@pde.PDEModel.delistMesh); 
    else      
        try
            [nodes,tet,cas,fas,eas,vas,Hmax,Hmin,Hgrad] = genmeshinternal(geom,Hmax,Hmin,Hgrad,Hface,Hedge,Hvertex,geomOrder);
        catch ex
            throwAsCaller(ex);
        end       
        assoc = pde.FEMeshAssociation(tet, cas, fas, eas, vas);     
        self.Mesh = pde.FEMesh(nodes,tet,Hmax,Hmin, Hgrad, geomOrder,assoc);
        msh = self.Mesh;
        addlistener(msh,'ObjectBeingDestroyed',@pde.PDEModel.delistMesh); 
    end
    
    function sizes = getValidLocalSizes(hc,num,err_index,err_format)
        sizes = zeros(1,num);
        if isempty(hc)
            return
        end
        if ~isvector(hc)
            error(message(err_format));
        end
        for i = 1:2:length(hc)
            id = hc{i};
            if ~isvector(id) || ~isreal(id) || ischar(id) || issparse(id(:)) || any(~isfinite(id(:))) ||...
                any(id < 1) || any(id > num)
                error(message(err_index));
            end
            sz = hc{i+1};
            if ~isreal(sz) || ~isscalar(sz) || ischar(sz) || sz < 0 || issparse(sz) || ~isfinite(sz)
                error(message('pde:pdeModel:invalidHlocal'));
            end
            sizes(id) = sz;
        end
    end
end
