function gm = importGeometry(filename, options)
% importGeometry - Creates a geometry object from the specified STL or STEP geometry file.
%    G = importGeometry(FILENAME) Constructs a geometry object from STL
%    or STEP data imported from a file. The STL data can represent a
%    3-D solid or a planar geometry. Planar geometry is supported by
%    STL format only. FILENAME is a character vector or string
%    specifying the name of the STL or STEP file ending with .stl, .STL,
%    .step, or .STEP and the path to it.
%    You can omit the path name if the file is in the
%    current directory or in a directory on the MATLAB path.
%    G is the new geometry object returned by importGeometry.
%
%    G = importGeometry(FILENAME,'MaxRelativeDeviation',d) lets you
%    control the accuracy of the geometry import from a STEP file by
%    specifying the relative sag. This parameter is ignored for STL files.
%
%    Example 1:  Import a planar geometry from an STL file. Plot the geometry.
%      gm = importGeometry('PlateHolePlanar.stl')
%      pdegplot(gm, 'EdgeLabels','on')
%
%    Example 2:  Import a 3-D geometry from an STL file. Plot the geometry.
%      gm = importGeometry('Block.stl')
%      pdegplot(gm, 'FaceLabels','on')
%
%    Example 3:  Import a 3-D geometry from a STEP file. Plot the geometry.
%      gm = importGeometry('BlockWithHole.step')
%      pdegplot(gm, 'FaceLabels','on')
%
%    See also pdegplot, pde.DiscreteGeometry, pde.DiscreteGeometry/addVertex

% Copyright 2020-2023 The MathWorks, Inc.

arguments
    filename {mustBeValidFileName}
    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
MaxRelativeDeviationDefault = 1;

if nargin ~= 1
    error(message('pde:pdeGeometryClass:InvalidNumInputArgs'));
end

% File exists and may be on the MATLAB path.
% Get full pathname and filename.
fid = fopen(filename, 'r');
filename = fopen(fid);
fclose(fid);

% Verify the file extension
[~,~,ext] = fileparts(filename);
if ~strcmpi(ext,'.stl') && ~strcmpi(ext,'.step') && ~strcmpi(ext,'.stp')
    error(message('pde:pdeGeometryClass:InvalidFileExtension'));
end

% apply deviation depending on the extension
if strcmpi(ext,'.stl')
    if isfield(options,'MaxRelativeDeviation')
        error(message('pde:pdeGeometryClass:STEPDeviationNotForStl'));
    end
    ST = matlab.internal.meshio.stlread(filename,4*eps('single'));
    if ST.ErrorCode > 0
        error(message('pde:pdeGeometryClass:InvalidFileFormat','STL'));
    end

    % check if the triangulation is valid and flip inconsistent normals
    try
        ST = checkAndFlipNormals(ST);
        checkMeshValidity(ST.Vertices',ST.Faces',[],ST.Normals);
    catch ME
        switch ME.identifier
            case 'pde:pdeGeometryClass:MeshHasGaps'
                error(message('pde:pdeGeometryClass:GeometryNotWatertight'));
            case  'pde:pdeGeometryClass:MeshIncompatible'
                error(message('pde:pdeGeometryClass:GeometryNonmanifold'));
            otherwise
                rethrow(ME)
        end
    end
    gm = geomFromMeshInternal(ST.Vertices',ST.Faces',[],ST.Normals,options.FeatureAngle);
else
    if isfield(options,'MaxRelativeDeviation')
        gm = pde.DiscreteGeometry(filename,options.MaxRelativeDeviation);
    else
        gm = pde.DiscreteGeometry(filename,MaxRelativeDeviationDefault);
    end
end
if ~options.AllowSelfIntersections
    checkForSelfIntersections(gm);
end
removeScratchEdges(gm);

end

function mustBeValidFileName(fn)
if ismissing(fn)
    error(message('pde:pdeGeometryClass:FileDoesNotExist',''));
end
fid = fopen(fn, 'r');
if (fid == -1)
    dirfile = dir(fn);
    if ~isempty(dirfile)
        if (length(dirfile) ~= 1)
            error(message('pde:pdeGeometryClass:FileIsDirectory'));
        else
            error(message('pde:pdeGeometryClass:FileOpenFailed'));
        end
    else
        error(message('pde:pdeGeometryClass:FileDoesNotExist',fn));
    end
else
    fclose(fid);
end
end


% check and flip inconsistent normals by swapping vertices of triangles
%
% Note 1: This makes triangle orientation consistent with the normals given
% in the STL file. It does not check that adjacent triangles have
% consistent orientation.
%
% Note 2: Orientation is not flipped if the original triangulation is
% watertight and the flipped one is not. It means that the imported normals
% are self-inconsistent while the triangulation is.
function ST = checkAndFlipNormals(ST)
if isempty(ST.Faces)
    error(message('pde:pdeGeometryClass:InvalidGeometryData'));
end
tri = triangulation(ST.Faces, ST.Vertices);
normals = faceNormal(tri);
product = dot(normals,ST.Normals,2);
tf = (product < -0.1);
if any(tf)
    faces = ST.Faces;
    faces(tf, [2 3]) = faces(tf, [3 2]);
    % flip only if it does not break watertight check
    if ~checkMeshWatertight(ST.Faces) || checkMeshWatertight(faces)
        ST.Faces = faces;
    end
end
end
