function [f] = fvector(node,elem,g0,ki)
%[f] = fvector(node,elem,g0,ki)
%
%    f: The force vector
% node: Nodal positions
% elem: Element connectivities
%   g0: The surface tension
%   ki: Used for anisotropy


%% fb the force on the two ends

% g0=30000;

fb=0.5;    % y direction

% gammaB=0.01;
% gammaS=0.2;

gammaB=0.1;
gammaS=0.5;

%Number of elements
nel=size(elem,1);

%Number of nodes
nnode=size(node,1);



   connec=zeros(nnode,3);
for ss=1:nnode
   index=1;
    for kk=1:nel
        if elem(kk,1)==ss || elem(kk,2)==ss
            connec(ss,index)=kk;
            index=index+1;
        end
    end
      
end
 
 

%Initialize the f vector
f=zeros(nel+3*nnode,1);


for i=1:nel;
    
    %The nodes on either side of the element
    node1=elem(i,1);
    node2=elem(i,2);
    
    %The nodal positions
    x1=node(node1,1);
    y1=node(node1,2);
    x2=node(node2,1);
    y2=node(node2,2);
    
    %Length of the element
    l=sqrt((y2-y1)^2+(x2-x1)^2);
    
    %The Sine and Cosine of the element angle
    S=(y2-y1)/l;
    C=(x2-x1)/l;
    
    %Find the forces at the node.
    %% added to select the right gamma
 if elem(i,3)==1;  %%the element is on boundary
       gamma=gammaB;
     
    elseif elem(i,3)==0; %% the element is on surface
        gamma=gammaS;
 
    else
        disp('error when selecting gamma');
 end
    
 %%
%     gamma=g0*(1+ki*abs(S));
    
    %Assemble the global f vector       
%         n1=elem(i,1);
%         n2=elem(i,2);
%         f(3*n1-2)=f(3*n1-2)+gamma*pi*(-l+(x1+x2)*C);
%         f(3*n1-1)=f(3*n1-1)+gamma*pi*((x1+x2)*S);
%         f(3*n2-2)=f(3*n2-2)+gamma*pi*(-l-(x1+x2)*C);
%         f(3*n2-1)=f(3*n2-1)+gamma*pi*(-(x1+x2)*S);
%         

   n1=elem(i,1);
        n2=elem(i,2);
        f(3*n1-2)=f(3*n1-2)+gamma*C;
        f(3*n1-1)=f(3*n1-1)+gamma*S;
        f(3*n2-2)=f(3*n2-2)+gamma*(-C);
        f(3*n2-1)=f(3*n2-1)+gamma*(-S);
        


    
end

%         f(3*1-1)=f(3*1-1)+0.8;
%         f(3*4-1)=f(3*4-1)+0.8;
%         
        
        
 for nn=1:nnode
   if nnz(connec(nn,:))==1 
      
        f(3*nn-1)=f(3*nn-1)+0.05;
   end
     
    
end

        
end

