function [mg, mesh] = Build_GMG_Data(mesh,Stiff, Mass, u_D, Levels)

% Build the geometric multigrid data on the mesh to produce the matrices
% corresponding to Stiff and Mass (stiffness and mass matrices)
% INPUT: mesh: the given initial mesh
%        Stiff: the matrix struct for producing the stiff matrix
%        Mass: the matrix struct for producing the mass matrix
%        u_D: the Dirichlet boundary condition
%        Levels: the levels for the multigrid method
% OUTPUT: mg: the multigrid data
%         mesh: the finest mesh
% Hehu Xie: 2015-09-12

%----initialize the MG data----
A = Matrix_Assemble(mesh,Stiff);
if isstruct(Mass)
    M = Matrix_Assemble(mesh,Mass);
    [A,M] = Poission_BoundaryTreatment(mesh,Stiff.L_Base,u_D,A,M);
    [A,M,bd] = Poission_BoundaryTreatment(mesh,Stiff.L_Base,u_D,A',M');
    mg = MGInit( A, M, bd);
else
    N = size(A,1);
    [A,tmp,bd] = Poission_BoundaryTreatment(mesh,Stiff.L_Base,u_D,A,sparse(N,1));
    mg = MGInit(A,[],bd);
end
%----deal with the boundary condition----------------------------------
Cur_Level = 2;
xv = sparse(size(A,1),1);
%==============Refine the mesh and build the multigrid system==============
while Cur_Level <= Levels
    Cur_Level
    [mesh,xv,Prolong] = refinement(mesh,xv,Stiff.L_Base);
    N(Cur_Level) = size(mesh.elem,1);
    %------First do the eigenvalue problem solving in the coarse space-----
    %-----build the stiff matrix ------------------------------------------
    A = Matrix_Assemble(mesh,Stiff);
    if isstruct(Mass)
        M = Matrix_Assemble(mesh,Mass);
        [A,M] = Poission_BoundaryTreatment(mesh,Stiff.L_Base,u_D,A,M);
        [A,M,bd] = Poission_BoundaryTreatment(mesh,Stiff.L_Base,u_D,A',M');
        mg = MGAddLevel(mg,A,M,Prolong,bd);
    else
        N = size(A,1);
        [A,tmp,bd] = Poission_BoundaryTreatment(mesh,Stiff.L_Base,u_D,A,sparse(N,1));
        mg = MGAddLevel(mg,A,[],Prolong,bd);
    end  
    Cur_Level = Cur_Level + 1;
end