function [g, m] = geometryFromMesh(self, nodes, elements, regions)
% geometryFromMesh - Creates a geometric model from a triangulated mesh
%    geometryFromMesh Creates a geometric model from a planar triangulation,
%    a volume triangulation or a surface triangulation that bounds a closed
%    volume. For planar and volume triangulated meshes, the mesh is retained
%    by default and can be deleted/regenerated as necessary.
%
%    geometryFromMesh(PDEM, Nodes, Elements) constructs a geometric model from
%    a mesh defined by Nodes and Elements. The variable Nodes is a matrix of
%    size 3-by-nnodes, where nnodes is the number of nodes in the mesh. The
%    variable Elements is a matrix of size mnpe-by-nelem where nelem is the
%    number of elements in the mesh and mnpe is the number of nodes per
%    element; 3 for a planar or surface mesh, 4 or 10 for a volume mesh. The
%    geometry created is assigned to the Geometry property of this class.
%    If the input is a planar or volume mesh, it is likewise assigned to the
%    Mesh property.
%
%    geometryFromMesh(PDEM, Nodes, Elements, ElementIdToRegionId)
%    constructs a multidomain geometric model from a mesh defined by Nodes
%    and Elements. Domain information is specified by ElementIdToRegionId,
%    a vector of length equal to the number of elements in the mesh.
%    ElementIdToRegionId(i) if the ID of the geometric region associated
%    with element i.
%
%    [G, MSH] = geometryFromMesh(__) returns handles to the geometry G and
%    mesh MSH. If the input mesh is not a volume mesh, MSH is empty, the mesh
%    can be generated using the generateMesh method.
%
%    Example 1: Construct a geometric block from a triangulated surface
%               The triangulation is obtained from the convex hull of a
%               meshgrid of points.
%      [x,y,z] = meshgrid(-2:4:2);
%      x = x(:); y = y(:); z = z(:);
%      K = convhull(x,y,z);
%      pdem = createpde();
%      geometryFromMesh(pdem, [x'; y'; z'], K')
%      pdegplot(pdem);
%
%    Example 2: Create a geometric model from a volume tetrahedral mesh
%      load tetmesh
%      % This loads variables X and tet into the workspace
%      % X represents the nodes and tet represents the elements
%      pdem = createpde();
%      geometryFromMesh(pdem, X', tet')
%      pdegplot(pdem);
%
%    Example 3: Create 2-D multidomain geometry from a mesh
%       load MotorMesh2D
%       %This loads variables nodes, elems, and eid2fid - a vector that
%       %returns the face ID for a given element ID index.
%       pdem = createpde;
%       pdem.geometryFromMesh(nodes,elems,eid2fid)
%       %Plot wireframe geometry
%       pdegplot(pdem, 'FaceLabels','on')
%
%       % Plot the mesh in a separate figure
%       % First find the element IDs associated with each face.
%       msh = pdem.Mesh;
%       eid1 = msh.findElements('region','face',1);
%       eid2 = msh.findElements('region','face',2);
%       eid3 = msh.findElements('region','face',3);
%       eid4 = msh.findElements('region','face',4);
%       eid5 = msh.findElements('region','face',5);
%       % Assign a different color to the elements in each subdomain.
%       figure('Color','white')
%       hold on
%       pdemesh(nodes,elems(:,eid3),'EdgeColor','magenta')
%       pdemesh(nodes,elems(:,eid5),'EdgeColor','blue')
%       pdemesh(nodes,elems(:,eid1),'EdgeColor','cyan')
%       pdemesh(nodes,elems(:,eid4),'EdgeColor','green')
%       pdemesh(nodes,elems(:,eid2),'EdgeColor','yellow')
%       axis off
%
%    Example 4: Create 3-D multidomain geometry from a mesh
%       load MultidomainMesh3D
%       %This loads variables nodes, elements, and ElementIdToRegionId
%       pdem = createpde();
%       geometryFromMesh(pdem, nodes, elements, ElementIdToRegionId);
%       pdegplot(pdem,'CellLabels','on');
%
%    Example 5: Create block with a cylindrical hole using an alphaShape.
%              Then create a geometry from the boundary triangulation
%              of the alphaShape.
%      Steps: Create a square grid of points and a radial grid of points
%              contained within it. Boolean these points using an alpha
%              shape, then sweep the result to 3D.
%      [xg, yg] = meshgrid(-3:0.25:3);
%      xg = xg(:);
%      yg = yg(:);
%      t = (pi/24:pi/24:2*pi)';
%      x = cos(t);
%      y = sin(t);
%      circShp = alphaShape(x,y,2);
%      in = inShape(circShp, xg,yg);
%      xg = [xg(~in); cos(t)];
%      yg = [yg(~in); sin(t)];
%      zg = ones(numel(xg),1);
%      % Sweep the points in Z
%      xg = repmat(xg,5,1);
%      yg = repmat(yg,5,1);
%      zg = zg*(0:.25:1);
%      zg = zg(:);
%      shp = alphaShape(xg,yg,zg, 0.5);
%      [tri, pts] = boundaryFacets(shp);
%      pdem = createpde();
%      geometryFromMesh(pdem,pts', tri');
%      pdegplot(pdem);
%      generateMesh(pdem);
%      figure
%      pdemesh(pdem)
%

%
%   See also pde.FEMesh, pde.DiscreteGeometry, pde.DiscreteGeometry/addVertex,
%            alphaShape, convhull, delaunay.

% Copyright 2015-2022 The MathWorks, Inc.

arguments
    self (1,1) pde.EquationModel {mustNotHaveGeometry}
    nodes (:,:) double {mustBeNumeric,mustBeNonempty,mustBeNonsparse,mustBeReal,mustBeFinite}
    elements (:,:) double {mustBeNumeric,mustBeNonempty,mustBeNonsparse,mustBeReal,mustBeFinite,mustBeInteger}
    regions double {mustBeNumeric,mustBeNonsparse,mustBeReal,mustBeFinite,mustBeInteger} = []
end

checkMeshValidity(nodes,elements,regions);

try
    if nargin < 4
        [g, m] = geomFromMeshInternal(nodes, elements);
    elseif isempty(regions)
        error(message('pde:pdeModel:InvalidSubdomainId'));
    else
        [g, m] = geomFromMeshInternal(nodes,elements,regions);
    end
catch ex
    throwAsCaller(ex);
end

self.Geometry = g;
self.Mesh = m;
if nargout == 2 && isempty(m)
    warning(message('pde:pdeModel:invalidDomainMesh'));
end

end

function mustNotHaveGeometry(pdem)
if ~isempty(pdem.Geometry)
    error(message('pde:pdeModel:noAssemSupport'));
end
end
