function [k,t]=pdeintrn(p,e,t)
%PDEINTRN Determine interior triangles in a triangular mesh.
%
%       Optionally update t with respect to subdomain number.

%       L. Langemyr 11-25-94.
%       Copyright 1994-2023 The MathWorks, Inc.
if nargout==1
    np = size(p,2);
    i  = find(e(6,:));              % turning left from boundary
    E  = sparse(e(1,i),e(2,i),1,np,np);
    i  = find(e(7,:));              % turning right from boundary
    E  = E | sparse(e(2,i),e(1,i),1,np,np);
    T  = sparse(t([1 2 3],:),t([2 3 1],:),ones(3,1)*(1:size(t,2)),np,np);
    [~,~,k] = find(T.*E);
    k1 = [];
    k  = unique(k);
    while length(k1)~=length(k)
        T1 = sparse(t([1 2 3],k),t([2 3 1],k),1,np,np)-E>0;
        [~,~,k2] = find(T1'.*T);
        k1 = k;
        k  = unique([k; k2]);
    end
else
    np = size(p,2);
    i  = find(e(6,:));              % turning left from boundary
    E  = sparse(e(1,i),e(2,i),e(6,i),np,np);
    i  = find(e(7,:));              % turning right from boundary
    E  = E + sparse(e(2,i),e(1,i),e(7,i),np,np);
    T  = sparse(t([1 2 3],:),t([2 3 1],:),ones(3,1)*(1:size(t,2)),np,np);
    [~,~,k] = find(T.*sign(E));
    [~,~,s] = find(sign(T).*E);
    k1 = [];
    [k,i,~] = unique(k);
    s = s(i);
    while length(k1)~=length(k)
        T1 = sparse(t([1 2 3],k),t([2 3 1],k),1,np,np)-sign(E)>0;
        T1 = T1'.*T;
        [~,~,k2] = find(T1);
        S = sparse(t([1 2 3],k),t([2 3 1],k),ones(3,1)*s',np,np);
        [~,~,s2] = find(S'.*sign(T1));
        k1 = k;
        [k,i,~] = unique([k;k2]);
        s(end+1:end+numel(s2)) = s2;
        s = s(i);
    end
    t(4,:) = 0;
    t(4,k) = s';
end
end