function [b,gb,Hessian,coord]=produce_Lagrange_basis(K)

%  Generate the base function for standard Lagrange $P_k$ element  
%
% INPUT:  the order of the polynomials
%
% OUTPUT: b_1: the basis function for $P_k^+$
%         gb_1: the corresponding gradient
%        
% Hehu Xie 9-12-2009
%--------------------------------------------------------------------------
syms x y A base v bubble coord mon_K_1 mon_K_2

%--------build the polynomial bases for the $P_K^+$
v=[];
%-----standard $P_k$
for i=0:K
    for j=0:K-i
        v=[v,x^i*y^j];
    end
end

%---------get the point coordination of the standard interpolation points--
director=[-1 1; 0 -1;1 0];
start_point=[1 0;0 1;0 0];
coord=[];
coord=[0 0; 1 0;0 1];
for i=1:K-1
    coord=[coord;start_point+i/K*director];
end

if K>2
    for j=1:K-2
        for i=1:K-1-j
            coord=[coord;[0 0]+[i/K, j/K]];
        end
    end
end

%--------------------------------------------------------------------------
dim=size(v,2);
dim_1=size(coord,1);
%-------------standard interpolation part----------------------------------
for i=1:dim_1
    for j=1:dim
        A(i,j)=double(subs(subs(v(j),x,coord(i,1)),y,coord(i,2)));
    end
end


%----------get the basis functions-----------------------------------------
b=v*inv(A);
%----------get the gradient basis functions--------------------------------
gb=[diff(b,x,1);diff(b,y,1)];
%---------Get the Hessian matrix function----------------------------------
Hessian = [diff(b,x,2);diff(diff(b,x,1),y,1);diff(diff(b,y,1),x,1);diff(b,y,2)];