function [f,fe] = f_global(elem, node, LM, g0, gI, gII, v)
global N_step
% deltaT
% pause
n_node=size(node,1);
n_elem=size(elem,1);
[temp, temp2]=size(node);
DOF=temp2-1;
n_eqn=DOF*n_node;

% local force vector fe
fe=zeros(4,n_node);
% n_elem
% pause
for e=1:n_elem;
    node1=elem(e,2);
    node2=elem(e,3);

    x1=node(node1,2);
    y1=node(node1,3);
    x2=node(node2,2);
    y2=node(node2,3);

    %The length of the element, sin, cos
    [l,angle] = lenangle(x1,y1,x2,y2);
    S=sin(angle);
    C=cos(angle);
    %
    %     S=(y2-y1)/l;
    %     C=(x2-x1)/l;


    % Anisotropic surface tension and diff value
    Temp = 300 *(1-0.2*cos(6*angle));%change the value of coefficiet of angle
    dTemp = 300*0.2*6*sin(2*angle);%change delta T also 
    gamma=Temp/100;
    deltagamma=dTemp/100;

    %     N_step
    %    pause
    %if (v > 300)
        %if (angle > 0 && angle< 3.14*2)
       % if (y1 > 0 && y2 > 0)
          % if (abs(y1-y2)<0.01)
            %gamma = 0.02*g0;
           % deltagamma =0.0;
       % end
%end

  %   end
  %  end

    fe(1,e)=gamma*C+deltagamma*(-S)+l*(gI-gII)/2*(-S);
    fe(2,e)=gamma*S+deltagamma*C+l*(gI-gII)/2*C;
    fe(3,e)=gamma*(-C)+deltagamma*S+l*(gI-gII)/2*(-S);
    fe(4,e)=gamma*(-S)+deltagamma*(-C)+l*(gI-gII)/2*C;
end

% assemble global force vector f
f=zeros(n_eqn,1); % initialize

[r,c]=size(LM);

for j=1:c; %n_elem
    for i=1:r; % r=4
        P=LM(i,j);
        f(P)=f(P)+fe(i,j);
    end
end


