function [LAM,N,Iter,Energy] = BEC_FEM_Multigrid(Max_Level,mesh)

% Use the multigrid (multilevel correction method) finite element method to solve 
% the eigenvalue problem: 
% $-\Delta u + W u +\zeta u^2 u =\lambda u $ 
% with the homogeneous boundary value condition
%  INPUT: Max_Level: the number of sequence of meshes which are produced by
%  the regular refinement from  the initial mesh
%         mesh: the input mesh (the initial mesh)
%  OUTPUT: LAM: the eigenvalue approximations,
%          N: the number of elements of the sequence of meshes
%          Iter: the nonlinear iteration times in each level of mesh
%          Energy: the energy in each level mesh
% Here, we use the simple self-consistent field ("fix-point") iteration
% without damping technique. For large $\ztea$, we suggest to use the
% damping technique to guarantee the convergence.
%
% Hehu Xie: 2015-02-29



if nargin < 2
    uniformrefine = 3;
    %--------unit square: triangle-----------------------------------------
    node = [0 0; 1 0; 1 1; 0 1; 0.5 0.5];
    elem = [ 5 1 2;
             5 2 3; 5 3 4; 5 4 1];
    Dirichlet = [1 2;2 3; 3 4; 4 1];
    Neumann=[];
    %-----Create the mesh--------------------------------------------------
    mesh=Initial_FEM_matlab(node,elem,Dirichlet,Neumann);
    for i=1:uniformrefine
        mesh = refinement(mesh);
    end
end
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
coord3 = quad_triangle(13);
order = 1;
Base = Lagrange_basis_function(order);
Stiff = Init_Matrix(1,1,Base,Base,[1 1 1], [1 1 1],@Discrete_Form_Stiff,coord3);
%--------------------------------------------------------------------------
Mass = Init_Matrix(1,1,Base,Base,[1 0 0], [1 0 0],@Discrete_Form_Mass,coord3);
%------------- Assemble matrices ------------------------------------------
A1 = Matrix_Assemble(mesh,Stiff);
M1 = Matrix_Assemble(mesh,Mass);
[A1,M1] = Poission_BoundaryTreatment(mesh,Base,@u_D,A1,M1);
%-----solve a linear eigenvalue to get the initial solution----------------
%---------do the nonlinear iteration---------------------------------------
max_It = 200;
Accuracy = 1e-6;
Low = 16;
Up = 60;
%--------------------------------------------------------------------------
[xv,lm,iresult] = sptarn(A1,M1,Low,Up);
xv = xv(:,1);
[v,ind] = max(abs(xv));
xv = xv/sign(xv(ind));
xv = xv/sqrt(xv'*M1*xv);
%----------to do the nonlinear iteration-----------------------------------
Cur_Level = 1;
Stiff = Init_Matrix(1,1,Base,Base,[1 1 1], [1 1 1],@Discrete_Form_Nonlinear_Stiff,coord3);
for kk = 1 : max_It      
    Stiff.Aux_FEM = Init_Aux_FEM(1,xv,Base,[1 0 0]);   
    xv_old = xv;
    A1 = Matrix_Assemble(mesh,Stiff);
    [A1,M1] = Poission_BoundaryTreatment(mesh,Base,@u_D,A1,M1);
    [xv,lm,iresult] = sptarn(A1,M1,Low,Up);
    lm = lm(1)
    xv = xv(:,1);
    [v,ind] = max(abs(xv));
    sgn = sign(xv(ind));
    xv = xv/sgn;
    xv = xv/sqrt(xv'*M1*xv);
    res = min(norm(xv-xv_old), norm(xv+xv_old))
    if res < Accuracy
        fprintf('Nonlinear iteration time: %d\n',kk);
        break;
    end   
end
Iter(:,Cur_Level) = [kk];
%------Compute the energy--------------------------------
Energy(Cur_Level,1) = energy(mesh,Base,xv,coord3,@coeff_varphi);
LAM(Cur_Level,1) = lm; 
N(1) = size(mesh.elem,1);
%---------Start the multilevel correction process--------------------------
%-----------Save the initial mesh------------------------------------------
mesh = Init_Multi_Mesh(mesh);
meshs = [mesh,mesh];
%--------------------------------------------------------------------------
%---------Multigrid eigenvalue solution------------------------------------
Num_Eig = 1;
while Cur_Level < Max_Level
    Cur_Level = Cur_Level + 1
    %------------propolongation---------------------
    [meshs(1),xv] = refinement(meshs(1),xv,Base);
    N(Cur_Level,1) = size(meshs(1).elem,1);
    M = Matrix_Assemble(meshs,Mass);
    %xv = xv/sqrt(xv'*M*xv);
    %----------------do the multilevel correction--------------------------
    Stiffhh = Init_Matrix(1,1,Base,Base,[1 1 1],[1 1 1],@Discrete_Form_Nonlinear_Stiff,coord3);    
    StiffHH = Init_Matrix(2,2,Base,Base,[1 1 1],[1 1 1],@Discrete_Form_Nonlinear_Stiff,coord3);    
    StiffHh = Init_Matrix(1,2,Base,Base,[1 1 1],[1 1 1],@Discrete_Form_Nonlinear_Stiff,coord3);    
    MassHh  = Init_Matrix(1,2,Base,Base,[1 0 0],[1 0 0],@Discrete_Form_Mass,coord3);
    MHh = Matrix_Assemble(meshs,MassHh);
    Stiffhh.Aux_FEM = Init_Aux_FEM(1,xv,Base,[1 0 0]);
    A = Matrix_Assemble(meshs,Stiffhh); 
    RHS = lm*M*xv;
    [A,RHS] = Poission_BoundaryTreatment(meshs(1),Base,@u_D,A,RHS);
    xv = A\RHS;
    xv = xv/sqrt(xv'*M*xv);
    %--------------------------------------------------------------------------
    for kk = 1 : max_It        
        StiffHH.Aux_FEM = Init_Aux_FEM(1,xv,Base,[1 0 0]);
        StiffHh.Aux_FEM = Init_Aux_FEM(1,xv,Base,[1 0 0]);
        Stiffhh.Aux_FEM = Init_Aux_FEM(1,xv,Base,[1 0 0]);
        AHH = Matrix_Assemble(meshs,StiffHH);
        AHh = Matrix_Assemble(meshs,StiffHh);
        A = Matrix_Assemble(meshs,Stiffhh);        
        bb = AHh*xv;        
        AA = [AHH,bb; bb',xv'*A*xv];
        bb = MHh*xv;
        MM = [M1,bb;bb',xv'*M*xv];
        [AA,MM] = Poission_BoundaryTreatment(meshs(2),Base,@u_D,AA,MM);
        [AA,MM] = Poission_BoundaryTreatment(meshs(2),Base,@u_D,AA',MM');
        %-----Eigenvalue solving-----------------
        [xv2,lm,iresult] = sptarn(AA,MM,Low,Up);        
        % get the new solution 
        lm = lm(1)
        xv2 = xv2(:,1);
        xv = Get_Solution_Multilevel(Base,meshs(2),xv2(1:end-Num_Eig,:),...
                                     Base,meshs(1),xv2(end-Num_Eig+1:end,:),xv); 
        %---------unitlization--------------------------------------------
        [v,ind] = max(abs(xv));
        sgn = sign(xv(ind));
        xv = xv/sgn;
        xv = xv/sqrt(xv'*M*xv);
        if min(norm(xv-StiffHH.Aux_FEM.Vec),norm(xv+StiffHH.Aux_FEM.Vec)) < Accuracy
            fprintf('Nonlinear iteration time: %d\n',kk);
            break;
        end        
    end
    Iter(Cur_Level,1) = kk;
    LAM(Cur_Level,1) = lm;                                  
   %------Compute the energy--------------------------------
   Energy(Cur_Level,1) = energy(mesh(1),Base,xv,coord3,@coeff_varphi);%energy_new(mesh,xv);
end
if N(end) < 2e5
    figure(1)
    showFeFunction(meshs(1),xv)
end
LAM
N
Iter
Energy
% print the convergence order of the eigenvalue approximations
%---use the extrapolation method to get the "accurate" eigenvalue
LAM_Exact = (4^order*LAM(end)-LAM(end-1))/(4^order-1);
LAM_Err = abs(LAM-LAM_Exact);
figure(2)
loglog(N,LAM_Err,'-o', N,0.75*LAM_Err(2)*N(2)^(order)*N.^(-order),'-.',...
       'LineWidth',3,'markersize',14)
word = ['slope=-', num2str(order)];
legend('eigenvalue errors',word);
set(gca, 'FontSize', 14)
title('Eigenvalue errors of the finite element method','fontsize',18)
xlabel('Number of elements','fontsize',18)
ylabel('Errors','fontsize',18)


%--------------------------------------------------------------------------
%  Given conditions
%--------------------------------------------------------------------------
%-----------------Sub functions-----------------------
function v = Discrete_Form_Stiff(VL,VR,X,Aux_FEM,Matrix)
v = (X(:,1).^2+X(:,2).^2).*VL(:,1).*VR(:,1) + VL(:,2).*VR(:,2) + VL(:,3).*VR(:,3);
%----the discrete form for the nonlinear problem-------------
function v = Discrete_Form_Nonlinear_Stiff(VL,VR,X,Aux_FEM,Matrix)
zeta = 1;
v = (zeta*Aux_FEM(:,1).^2 +X(:,1).^2+X(:,2).^2).*VL(:,1).*VR(:,1) + VL(:,2).*VR(:,2) + VL(:,3).*VR(:,3);
%function v = Discrete_Form_Nonlinear_Tensor(VL,VR,V3,X,Aux_FEM,Matrix)
%zeta = 50;
%v = (zeta*Aux_FEM(:,1).^2+X(:,1).^2+X(:,2).^2).*VL(:,1).*VR(:,1).*V3(:,1);
%-----------------Mass Matrix---------------------
function v = Discrete_Form_Mass(VL,VR,X,Aux_FEM,Matrix)
v = VL(:,1).*VR(:,1);
%-----------Coefficient function degree 0----------------------------------
function v=coeff_varphi(x)
%--------Model problem-----------------------------------------------------
%-------Example 1----------------------------------------------------
v = x(:,1).^2+x(:,2).^2;
%--------------------------------------------------------------------------
function v=u_D(x)
v = sparse(size(x,1),1);
%---------------------------------------------------------------------------
function v=u_N(x)
v = sparse(size(x,1),1);