function checkMeshValidity(nodes,elements,regions,normals)
% Internal function to check whether the input data represents a valid mesh
% Copyright 2021-2023 The MathWorks, Inc.

% check node and element dimensions
if size(nodes,1) < 2 || size(nodes,1) > 3
    error(message('pde:pdeModel:NodesNot2dOr3d'));
end
if isempty(elements)
    error(message('pde:pdeGeometryClass:InvalidGeometryData'));
end

% check node ids
if min(elements(:)) < 1 || max(elements(:)) > size(nodes,2)
    error(message('pde:pdeModel:NodeIdOutOfRange'));
end

% classify triangulation
is_surface = size(elements,1) == 3 || size(elements,1) == 6;
is_volume = size(elements,1) == 4 || size(elements,1) == 10;
is_planar = size(nodes,1) == 2;
if (nargin > 3)
    is_planar = all(max(normals,[],1)-min(normals,[],1) < 2*eps('single'));
end
if (~is_surface && ~is_volume) || (is_volume && is_planar)
    error(message('pde:pdeModel:UnsupportedElementTopology'));
end
is_multiregion = ~isempty(regions);

% validate region data
if is_multiregion
    if (is_volume || is_planar)
        if ~isvector(regions) || numel(regions) ~= size(elements,2)
            error(message('pde:pdeModel:InvalidSubdomainId'));
        end

        % must be compact set of positive integers
        urs = unique(regions);
        if min(urs) < 1 || max(urs) > numel(urs)
            error(message('pde:pdeModel:InvalidSubdomainId'));
        end
    else
        error(message('pde:pdeModel:InvalidSubdomainId'));
    end
end

% check validity of single-domain mesh connectivity
if ~is_multiregion
    % check whether 3-D surface mesh is watertight
    if is_surface && ~is_planar && ~checkMeshWatertight(elements')
        error(message('pde:pdeGeometryClass:MeshHasGaps'));
    end
    % check whether triangle mesh is two-manifold
    if is_surface
        checkMeshIncompatible(elements(1:3,:)');
    end
    % check whether 1-domain solid mesh boundary is two-manifold
    if is_volume
        wos = warning('off','MATLAB:triangulation:PtsNotInTriWarnId');
        ocu = onCleanup(@()warning(wos));
        tr = triangulation(elements(1:4,:)',nodes');
        checkMeshIncompatible(freeBoundary(tr));
    end
end

% check validity of multi-domain mesh connectivity
if is_multiregion && is_volume
    % check whether each mesh domain boundary is two-manifold
    wos = warning('off','MATLAB:triangulation:PtsNotInTriWarnId');
    ocu = onCleanup(@()warning(wos));
    urs = unique(regions);
    for i = numel(urs)
        tr = triangulation(elements(1:4,regions == urs(i))',nodes');
        checkMeshIncompatible(freeBoundary(tr));
    end
end
end
