function I= fem_num_integral_loc(f,X,E,j)


% f----function handel to be integrated over
% X-----finite element nodes
% E----finite element conectivity matrix

    
    % get local nodes 
        nodes_local = E (:,j)';
        x_loc= X(nodes_local);
        a= x_loc (1);
        b= x_loc (2);
        
        % set quadrature nodes and weights 
        s = 1/ sqrt(3);
        X1= ((1+s)*a + (1-s)*b)/2; 
        X2= ((1-s)*a + (1+s)*b)/2;
        W1=0.5; W2= 0.5;
        
        I= (b-a)* (W1* f(X1) + W2* f(X2));  % update the integral
        
    
end