function F

global F vec_p vec_q norm_p norm_q Nml_Vct Area1 IEN LM ID B BNDY3 BC3

Sz1=size(IEN);

% compute the f (driving force induced by Gamma_s)
Gamma_s=1;
g=0;
Fs_3=zeros(B,1);Fs_1=zeros(B,1);Fs_2=zeros(B,1);
g3=g/3;

% element type 12 x 13  (p cross q)
for e=1:Sz1(1)
    vec_r(e,:)=vec_p(e,:)-vec_q(e,:);
    norm_r(e)=norm(vec_r(e,:),2);
   
    Nmpq_2=Area1(e,1);
    d1=-Gamma_s*(vec_p(e,1))/norm_p(e)^2;
    d2=-Gamma_s*(vec_p(e,2))/norm_p(e)^2;
    d3=-Gamma_s*(vec_p(e,3))/norm_p(e)^2;
    d4=-Gamma_s*(vec_q(e,1))/norm_q(e)^2;
    d5=-Gamma_s*(vec_q(e,2))/norm_q(e)^2;
    d6=-Gamma_s*(vec_q(e,3))/norm_q(e)^2;
    
    f(1,e)=Nmpq_2*(-d1-d4);
    f(2,e)=Nmpq_2*(-d2-d5);
    f(3,e)=Nmpq_2*(-d3-d6);
    f(4,e)=Nmpq_2*d1;
    f(5,e)=Nmpq_2*d2;
    f(6,e)=Nmpq_2*d3;
    f(7,e)=Nmpq_2*d4;
    f(8,e)=Nmpq_2*d5;
    f(9,e)=Nmpq_2*d6;
    
    for a=1:9
        P=LM(a,e);
        if P>0
           Fs_1(P,1)=f(a,e)+Fs_1(P,1);
       end
    end
end    


% element type 23 x 21  (-r cross -p)
for e=1:Sz1(1)

    Nmpq_2=Area1(e,1);
    d1=-Gamma_s*(-vec_r(e,1))/norm_r(e)^2;
    d2=-Gamma_s*(-vec_r(e,2))/norm_r(e)^2;
    d3=-Gamma_s*(-vec_r(e,3))/norm_r(e)^2;
    d4=-Gamma_s*(-vec_p(e,1))/norm_p(e)^2;
    d5=-Gamma_s*(-vec_p(e,2))/norm_p(e)^2;
    d6=-Gamma_s*(-vec_p(e,3))/norm_p(e)^2;
    
    f(1,e)=Nmpq_2*d4;
    f(2,e)=Nmpq_2*d5;
    f(3,e)=Nmpq_2*d6;
    f(4,e)=Nmpq_2*(-d1-d4);
    f(5,e)=Nmpq_2*(-d2-d5);
    f(6,e)=Nmpq_2*(-d3-d6);
    f(7,e)=Nmpq_2*d1;
    f(8,e)=Nmpq_2*d2;
    f(9,e)=Nmpq_2*d3;
    
    for a=1:9
        P=LM(a,e);
        if P>0
           Fs_2(P,1)=f(a,e)+Fs_2(P,1);
       end
    end
end

% element type 31 x 32  (-q cross r)
for e=1:Sz1(1)
    Nmpq_2=Area1(e,1);
    d1=-Gamma_s*(-vec_q(e,1))/norm_q(e)^2;
    d2=-Gamma_s*(-vec_q(e,2))/norm_q(e)^2;
    d3=-Gamma_s*(-vec_q(e,3))/norm_q(e)^2;
    d4=-Gamma_s*(vec_r(e,1))/norm_r(e)^2;
    d5=-Gamma_s*(vec_r(e,2))/norm_r(e)^2;
    d6=-Gamma_s*(vec_r(e,3))/norm_r(e)^2;
        
    f(1,e)=Nmpq_2*d1;
    f(2,e)=Nmpq_2*d2;
    f(3,e)=Nmpq_2*d3;
    f(4,e)=Nmpq_2*d4;
    f(5,e)=Nmpq_2*d5;
    f(6,e)=Nmpq_2*d6;
    f(7,e)=Nmpq_2*(-d1-d4);
    f(8,e)=Nmpq_2*(-d2-d5);
    f(9,e)=Nmpq_2*(-d3-d6);
    
    for a=1:9
        P=LM(a,e);
        if P>0
           Fs_3(P,1)=f(a,e)+Fs_3(P,1);
       end
    end
end

% phase change terms (free energy)
Fg=zeros(B,1);
for e=1:Sz1(1)
    g3l_=-g3*Nml_Vct(e,1);
    g3m_=-g3*Nml_Vct(e,2);
    g3n_=-g3*Nml_Vct(e,3);
    
    f(1,e)=Nmpq_2*g3l_;
    f(2,e)=Nmpq_2*g3m_;
    f(3,e)=Nmpq_2*g3n_;
    f(4,e)=Nmpq_2*g3l_;
    f(5,e)=Nmpq_2*g3m_;
    f(6,e)=Nmpq_2*g3n_;
    f(7,e)=Nmpq_2*g3l_;
    f(8,e)=Nmpq_2*g3m_;
    f(9,e)=Nmpq_2*g3n_;    

    for a=1:9
        P=LM(a,e);
        if P>0
           Fg(P,1)=f(a,e)+Fg(P,1);
       end
    end
end

Fs=1/3*(Fs_1+Fs_2+Fs_3)+Fg;

% compute Gamma_b
Gamma_b=1.41;
Sz2=size(BNDY3);
vec_z=[0 0 1]';

% add Gamma_b effect on boundary 3
Fb3=zeros(B,1);
IEN_f3(1,:)=BNDY3(1:Sz2(2)-1);
IEN_f3(2,:)=BNDY3(2:Sz2(2));
for eb3=1:Sz2(2)-1
    r3(1,eb3)=BC3(eb3+1,2)-BC3(eb3,2);
    r3(2,eb3)=BC3(eb3+1,3)-BC3(eb3,3);
    r3(3,eb3)=BC3(eb3+1,4)-BC3(eb3,4);
    norm_r3(eb3)=norm(r3(:,eb3),2);
    
    b3(:,eb3)=cross(cross(vec_z,r3(:,eb3)),r3(:,eb3));
    norm_b3(eb3)=norm(b3(:,eb3),2);
    Nml_b3(:,eb3)=b3(:,eb3)/norm_b3(eb3);
    
    LM_f3(1:3,eb3)=ID(:,IEN_f3(1,eb3));
    LM_f3(4:6,eb3)=ID(:,IEN_f3(2,eb3));
    
    CoefB3=1/2*Gamma_b*norm_r3(eb3);
    fb3(1:3,eb3)=CoefB3*Nml_b3(:,eb3);
    fb3(4:6,eb3)=CoefB3*Nml_b3(:,eb3);
    for a=1:6
        P=LM_f3(a,eb3);
        if P>0
            Fb3(P,1)=fb3(a,eb3)+Fb3(P,1);
        end
    end
end


 F=Fs+Fb3;

disp('F OK!')