function [mesh,ed4n]=Initial_FEM_matlab(node,elem,Dirichlet,Neumann)

% INITIAL_FEM_MATLAB generates initial mesh and buile its data structure.
% 
% USAGE
%    mesh = Initial_FEM_matlab(node,elem,Dirichlet,Neumann)
% INPUT 
%       node:  array of nodes (required) 
%       elem:  array of element indices (required)
%  Dirichlet:  Dirichlet boundary edges (required)
%    Neumann:  Neumann boundary edges (optional)
%    Dirichlet and Neumann stores the start node and end node for the
%    boundary edges
%
% OUTPUT
%    mesh:  current mesh
%    ed4n: the edge numbers for each node (for some special subroutines)
%
% Hehu Xie 05-17-2009

%--------------------------------------------------------------------------
% Generate opitional data
%--------------------------------------------------------------------------
tic;
fprintf('----------Generate the mesh data---------------------\n')
%---------------------- get the numbers of nodes and elements -------------          
nrElems = size(elem,1);
nrNodes = size(node,1);
%--------------------- get element number for eache node (e4n)--------------------
I = [elem(:,2);elem(:,3);elem(:,1)];
J = [elem(:,3);elem(:,1);elem(:,2)];
S = [1:nrElems;1:nrElems;1:nrElems];
e4n = sparse(I,J,S);
%--------------------- get edge number for each number (ed4n)-----------------------
ed4n = triu(e4n+e4n');
[I,J] = find(ed4n);
nrEdges = length(I);
ed4n = sparse(I,J,1:nrEdges,size(e4n,1),size(e4n,1));
ed4n = ed4n + ed4n';
[I,J,S] = find( triu(ed4n) );
n4ed( S,: ) = [I,J];
%-----------Initialize the boundary condition------------------------------
Nr_bd_Dirichlet = size(Dirichlet,1);
Dirichlet_edge = [];
if Nr_bd_Dirichlet > 0
    Dirichlet_edge = diag(ed4n(Dirichlet(:,1),Dirichlet(:,2)));
end
%--------------------------------------------------------------------------
Neumann_edge = [];
if nargin >= 4
    Nr_bd_Neumann = size(Neumann,1);
    if Nr_bd_Neumann > 0
        Neumann_edge = diag(ed4n(Neumann(:,1),Neumann(:,2)));
    end
else
    Neumann=[];       
end

%------------- get n4ed ---------------------------------------------------
[I,J,S] = find(triu(ed4n) );
n4ed( S,: ) = [I,J];
%----------get the edge information for each element-----------------------
%edge = [diag(ed4n(elem(:,2),elem(:,3))),diag(ed4n(elem(:,3),elem(:,1))),...
%diag(ed4n(elem(:,1),elem(:,2)))];
%----Use the following subroutine to save the memory cost------------------
edge = [rowaddr(ed4n,elem(:,[2 3 1]),elem(:,[3 1 2]))];
%type = ones(nrNodes,1);
%--------------------------------------------------------------------------
% Generate mesh data structure
%--------------------------------------------------------------------------
mesh = struct('node',node, ...
              'elem',elem(:,1:3), 'edge',edge,'node4edge',n4ed, ...
              'Dirichlet_node',Dirichlet,...
              'Dirichlet_edge',Dirichlet_edge,...
              'Neumann_node',Neumann,...
              'Neumann_edge',Neumann_edge,...
              'type',[]);   
toc               
%--------------------------------------------------------------------------
% End of function INITMESH
%--------------------------------------------------------------------------