function [E,T2E,Ed,En]=estendi_mesh(p,t,distD, distN)

% Funzione che estende la triangolazione t
% 
% INPUT:
% t = triangolazione 
% p = coordinate dei punti della triangolazione t
% distD = funzione distanza che restituisce 0 in corrispondenza del bordo
%         di Dirichlet
% distN = funzione distanza che restituisce 0 in corrispondenza del bordo
%         di Neumann
% 
% OUTPUT:
% E = matrice che contiene i lati orientati (ovvero gli indici dei vertici
%         del lato)
% T2E = matrice che alla riga i-esima contiene un'estensione della matrice
%         t: sono stati aggiunti gli indici dei lati opposti ai vertici e 
%         l'orientamento delle normali
% Ed = indici dei lati del bordo di Dirichlet
% En = indici dei lati del bordo di Neumann
% 
% 
% Versione del codice rivista da Giovanni Varetto (A.A. 2023-2024) a
% partire da un codice pre-esistente e contenente errori. In questa verione
% si usano le funzioni distD e distN per trovare (a meno di una tolleranza
% tol) i punti sul bordo ed i relativi lati.

T2E = t;
Ed = [];
En = [];
tol = 10 ^ -4;

for i=1:size(t,1)

    M=[t(i,2),t(i,3);
        t(i,3),t(i,1);
        t(i,1),t(i,2)];


    if i==1
        E(1:3,:) = M;
        T2E(i,3+[1:3])=[1:3];
        T2E(i,6+[1:3])=1;
        dim = 3;

        for  yy=1:3
            if abs(distD(p(E(yy,1),:))) + abs(distD(p(E(yy,2),:))) < tol
                Ed = [Ed yy];
            end

            if abs(distN(p(E(yy,1),:))) + abs(distN(p(E(yy,2),:))) < tol
                En = [En yy];
            end
        end
    else
        for n=1:3
            [ii]=find(E(:,1)==[M(n,2)]);
            flag=0;
            if ~isempty(ii)
                jj = find(E(ii,2)==M(n,1));
                if ~isempty(jj)
                    E(ii(jj),2)==M(n,1);
                    T2E(i,n+3)=ii(jj);
                    T2E(i,n+6) = -1;
                    flag =1;

                    yy=ii(jj);
                    if abs(distD(p(E(yy,1),:))) + abs(distD(p(E(yy,2),:))) < tol
                        Ed = [Ed yy];
                    end

                    if abs(distN(p(E(yy,1),:))) + abs(distN(p(E(yy,2),:))) < tol
                        En = [En yy];
                    end
                else
                    flag = 0;
                end
            end

            if flag ==0
                dim=dim+1;
                T2E(i,n+3) = dim;
                T2E(i,n+6) = 1;
                E(dim,:) = [M(n,1) M(n,2)];

                  yy=dim;
                    if abs(distD(p(E(yy,1),:))) + abs(distD(p(E(yy,2),:))) < tol
                        Ed = [Ed yy];
                    end

                    if abs(distN(p(E(yy,1),:))) + abs(distN(p(E(yy,2),:))) < tol
                        En = [En yy];
                    end

            end

        end

    end

Ed = unique(Ed);
En = unique(En);
end