function [Free_deg3,Dim]=DOF_Mangement(mesh,DOF3,Scale3)

% DOF_Mangement.m produce the degree of freedom corresponding to the degree
% freedom description DOF
% INPUT: MESH, the mesh used in the finite element method
%        DOF, the DOF description, DOF(1) denote the number of the degree 
%   of freedom  at each vertex, DOF(2), denote the number of the degree of
%   freedom on each edge and DOF(3) denote the number of the degree of
%   freedom in each element. For example for the P_3 element, the DOF is 
%  DOF=[1 2 1]
%  OUTPUT: FREE_DEG, the numbering for each degree of freedom
% ---------------Hehu Xie, 20-05-2009------------------------------------

%------set the initial values for output-----------------------------------
Free_deg3 = [];
Dim = 0;
if nargin == 2
    Scale3 = 0;
end


%------get the information of the mesh-------------------------------------
Nr_node=size(mesh.node,1);
Nr_edge=size(mesh.node4edge,1);
Nr_elem=size(mesh.elem,1);
%----get the number of triangles-------------------------------------------
% II3 = find(mesh.elem(:,4)<0);
% Nr_elem3 = length(II3);
% II4 = setdiff([1:Nr_elem]',II3);
% Nr_elem4 = length(II4);
II3 = 1 : Nr_elem;
%----------compute the freedom degree vector-------------------------------
Dim = Nr_node*DOF3(1)+Nr_edge*DOF3(2)+Nr_elem*DOF3(3);
if Nr_elem > 0
    Free_base = [kron(Nr_node*[0:DOF3(1)-1],ones(1,3)),...
        kron(Nr_node*DOF3(1)+Nr_edge*[0:DOF3(2)-1],ones(1,3)),...
        kron(Nr_node*DOF3(1)+Nr_edge*DOF3(2)+Nr_elem*[0:DOF3(3)-1],ones(1,1))];
    
    Free_deg3 = ones(Nr_elem,1)*Free_base+...
                [kron(ones(1,DOF3(1)),mesh.elem(:,1:3)),...
                 kron(ones(1,DOF3(2)),mesh.edge(:,1:3)),...
                 kron(ones(1,DOF3(3)),[1:Nr_elem]')];
     %----adjust the egde DOF with the edge direction----------------------
    if (DOF3(2) > 1)&(Scale3 >= 0)
        node4edge = [2,3,1];
        invDOF = DOF3(2)-1:-1:0;
        for i=1:3
            Index       = find(mesh.elem(:,node4edge(i)) ~= mesh.node4edge(mesh.edge(:,i),1));        
            Free_deg3(Index,3*DOF3(1)+ i+3*[0:DOF3(2)-1]) = Free_deg3(Index,3*DOF3(1)+ i+3*invDOF);
        end
    end 
end