function [E,E3,E4] = energy(mesh,Base,solu,coord3,coeff_varphi)
% Compute the energy for the BEC problem: 
% -\Delta u + Wu +\zeta u^2u = \lambda u 
% for the given solution solu with the basis function Base. 
%
% INPUT: mesh, the finite element mesh
%        Base: the basis of the finite element space
%        solu: the given finite element solution
%        coord3: the quadrature information
%        coeff_varphi: the funtion \zeta for the BEC problem
%
% OUTPUT: ground state energy
%
%


%--------get the element matrix size---------------------------------------
Nr_int3=size(coord3,1);
%--------compute the base value at the integral points---------------------
Nr_elem=size(mesh.elem,1);
%----------compute the freedom degree vector-------------------------------
[Free_deg1_3,Dim1] = DOF_Mangement(mesh,Base.DOF_3);
%----compute the information of the mesh ----------------------------------
[area,inv_transf]=Transform_information(mesh,1,1);
%------compute the first item(gradient)----------------
E1 = 0;
m1 = size(Base.gb_3([0,0]),2);
for kk = 1 : Nr_int3
    B_1 = Base.gb_3([0 0]*coord3(kk,1)+[1 0]*coord3(kk,2)+[0 1]*(1-coord3(kk,1)-coord3(kk,2)));
    %---------compute the gradient of the approximation-------------
    v = zeros(Nr_elem,2);
    for i = 1 : m1
        v = v + [(inv_transf(:,1,1)*B_1(1,i)+inv_transf(:,2,1)*B_1(2,i)).*solu(Free_deg1_3(:,i)),...
            (inv_transf(:,1,2)*B_1(1,i)+inv_transf(:,2,2)*B_1(2,i)).*solu(Free_deg1_3(:,i))];
    end
    E1 = E1 + coord3(kk,3)*sum(sum(v.^2,2).*area);
end
%E1 = 0.5*E1;

%------compute the second item(W)---------------
E2 = 0;
%[area]=Transform_information(mesh,1,0,3);
for kk = 1 : Nr_int3
    coeff_matrix=coeff_varphi((1-coord3(kk,1)-coord3(kk,2))*mesh.node(mesh.elem(:,1),1:2)+...
                       coord3(kk,1)*mesh.node(mesh.elem(:,2),1:2)+...
                       coord3(kk,2)*mesh.node(mesh.elem(:,3),1:2));
    B_1 = Base.b_3([0 0]*coord3(kk,1)+[1 0]*coord3(kk,2)+[0 1]*(1-coord3(kk,1)-coord3(kk,2)));
    %----------compute the error-----------------------------------
    E2 = E2 + coord3(kk,3)*sum((coeff_matrix.*((solu(Free_deg1_3)*B_1').^2)).*area);
end

%------compute the third item(nonlinear)---------------
E3 = 0;
%[area]=Transform_information(mesh,1,0,3);
for kk = 1: Nr_int3
    B_1 = Base.b_3([0 0]*coord3(kk,1)+[1 0]*coord3(kk,2)+[0 1]*(1-coord3(kk,1)-coord3(kk,2)));
    %----------compute the error-----------------------------------
    E3 = E3 + coord3(kk,3)*sum((solu(Free_deg1_3)*B_1').^4.*area);
end
E3 = 0.5*E3;

%------compute the norm item---------------
E4 = 0;
%[area]=Transform_information(mesh,1,0,3);
for kk = 1: Nr_int3

    B_1 = Base.b_3([0 0]*coord3(kk,1)+[1 0]*coord3(kk,2)+[0 1]*(1-coord3(kk,1)-coord3(kk,2)));
    %----------compute the error-----------------------------------
    E4 = E4 + coord3(kk,3)*sum(((solu(Free_deg1_3)*B_1').^2).*area);
end

E=(E1+E2+E3)/E4;

% 
% function v=coeff_varphi(x)
% %--------Model problem-----------------------------------------------------
% %v = zeros(size(x,1),1);
% %-------Example 1----------------------------------------------------
% v = x(:,1).^2+x(:,2).^2;