function [e4edge,Local_Num,Edge3,Edge4] = element4edge(mesh,bd_edge)

%  find the element number for the edge and to produce the information of 
%  element4edge
% INPUT: mesh: the computing finite element mesh
% OUT:  element4edge, the global element number for each edge
%        Local_Num: the local edge number in the corresponding element
%
% HEHU XIE, 2010-12-27

%---------------------- get the numbers of nodes and elements -------------          
nrElems = size(mesh.elem,1);
nrNodes = size(mesh.node,1);
%------get the number of the triangles-------------------------------------
II3 = find(mesh.elem(:,4)<0);
nrElems3 = length(II3);
% ------------treat the quadrilateral case --------------------------------
II4 = setdiff([1:nrElems]',II3);
nrElems4 = length(II4);
%--------------------- get element4edge------------------------------------
if nrElems4 > 0
    I = [II3,II3,II3,II4,II4,II4,II4];
    J = [mesh.edge(II3,1:3),mesh.edge(II4,1:4)];
    S = [II3,II3,II3,II4,II4,II4,II4];
else
    I = [II3,II3,II3];
    J = [mesh.edge(II3,1:3)];
    S = [II3,II3,II3];
end

e4ed = sparse(I,J,S);
Edge3=[];
Edge4=[];
%--------produce the e4ed--------------------------------------------------
num_edge = size(bd_edge,1);
e4edge = sparse(num_edge,1);
for k = 1:num_edge
    index = find(e4ed(:,bd_edge(k)));
    e4edge(k,:) = index;
    if mesh.elem(index,4) < 0
        Edge3 = [Edge3;k];
    else
        Edge4 = [Edge4;k];
    end
    Local_Num(k,:) = find(~(mesh.edge(index,:)-bd_edge(k)*ones(1,4)));
end