function H

global H vec_p vec_q norm_p norm_q Nml_Vct Area1 LM IEN B Node

% compute the vectors and normal vectors of the element surfaces
% vector p is "node 1 to node 2", vector q is "node 1 to node 3", 
% Nml_Vct is the normal unit vector perpendicular to the element surface
Sz1=size(IEN);   % number of elements
for e=1:Sz1(1)
    vec_p(e,1)=Node(IEN(e,3),2)-Node(IEN(e,2),2);
    vec_p(e,2)=Node(IEN(e,3),3)-Node(IEN(e,2),3);
    vec_p(e,3)=Node(IEN(e,3),4)-Node(IEN(e,2),4);
    vec_q(e,1)=Node(IEN(e,4),2)-Node(IEN(e,2),2);
    vec_q(e,2)=Node(IEN(e,4),3)-Node(IEN(e,2),3);
    vec_q(e,3)=Node(IEN(e,4),4)-Node(IEN(e,2),4);
    CRS1(e,:)=cross(vec_p(e,:),vec_q(e,:));
    Area1(e,1)=1/2*norm(CRS1(e,:),2);
    Nml_Vct(e,:)=CRS1(e,:)/2/Area1(e,1);
    norm_p(e)=norm(vec_p(e,:),2);
    norm_q(e)=norm(vec_q(e,:),2);
    
    Nml_p(e,:)=vec_p(e,:)/norm_p(e);
    
    X1(e)=0; Y1(e)=0;
    X2(e)=norm_p(e); Y2(e,2)=0;
    X3(e)=dot(vec_q(e,:),Nml_p(e,:)); Y3(e)=norm((vec_q(e,:)-dot(vec_q(e,:),Nml_p(e,:))*Nml_p(e,:)),2);
    
    J(e)=(X2(e)-X1(e))*(Y3(e)-Y1(e))-(X3(e)-X1(e))*(Y2(e)-Y1(e));
    
end
% vec_p
% vec_q
% CRS1
% Area1
% Nml_Vct
% norm_p
% norm_q

% build the global H matrix
H=zeros(B,B);
L_mo=10^2;     % mobility

for e=1:Sz1(1,1)
    
    % H_NN is the local H matrix, the following is the integrated form
    H_NN(1,1,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,1)*J(e);
    H_NN(1,2,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(1,3,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    H_NN(1,4,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,1)*J(e);
    H_NN(1,5,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(1,6,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    H_NN(1,7,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,1)*J(e);
    H_NN(1,8,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(1,9,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    
    H_NN(2,1,e)=H_NN(1,2,e);
    H_NN(2,2,e)=1/L_mo*1/12*Nml_Vct(e,2)*Nml_Vct(e,2)*J(e);
    H_NN(2,3,e)=1/L_mo*1/12*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e);
    H_NN(2,4,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(2,5,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,2)*J(e);
    H_NN(2,6,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e);
    H_NN(2,7,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(2,8,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,2)*J(e);
    H_NN(2,9,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e);
    
    H_NN(3,1,e)=H_NN(1,3,e);
    H_NN(3,2,e)=H_NN(2,3,e);
    H_NN(3,3,e)=1/L_mo*1/12*Nml_Vct(e,3)*Nml_Vct(e,3)*J(e);
    H_NN(3,4,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    H_NN(3,5,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e);
    H_NN(3,6,e)=1/L_mo*1/24*Nml_Vct(e,3)*Nml_Vct(e,3)*J(e);
    H_NN(3,7,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    H_NN(3,8,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e);
    H_NN(3,9,e)=1/L_mo*1/24*Nml_Vct(e,3)*Nml_Vct(e,3)*J(e);  
    
    H_NN(4,1,e)=H_NN(1,4,e);
    H_NN(4,2,e)=H_NN(2,4,e);
    H_NN(4,3,e)=H_NN(3,4,e);
    H_NN(4,4,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,1)*J(e);
    H_NN(4,5,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(4,6,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    H_NN(4,7,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,1)*J(e);
    H_NN(4,8,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(4,9,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    
    H_NN(5,1,e)=H_NN(1,5,e);
    H_NN(5,2,e)=H_NN(2,5,e);
    H_NN(5,3,e)=H_NN(3,5,e);
    H_NN(5,4,e)=H_NN(4,5,e);
    H_NN(5,5,e)=1/L_mo*1/12*Nml_Vct(e,2)*Nml_Vct(e,2)*J(e);
    H_NN(5,6,e)=1/L_mo*1/12*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e); 
    H_NN(5,7,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(5,8,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,2)*J(e);
    H_NN(5,9,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e);
    
    H_NN(6,1,e)=H_NN(1,6,e);
    H_NN(6,2,e)=H_NN(2,6,e);
    H_NN(6,3,e)=H_NN(3,6,e);
    H_NN(6,4,e)=H_NN(4,6,e);
    H_NN(6,5,e)=H_NN(5,6,e);
    H_NN(6,6,e)=1/L_mo*1/12*Nml_Vct(e,3)*Nml_Vct(e,3)*J(e);
    H_NN(6,7,e)=1/L_mo*1/24*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    H_NN(6,8,e)=1/L_mo*1/24*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e);
    H_NN(6,9,e)=1/L_mo*1/24*Nml_Vct(e,3)*Nml_Vct(e,3)*J(e);
    
    H_NN(7,1,e)=H_NN(1,7,e);
    H_NN(7,2,e)=H_NN(2,7,e);
    H_NN(7,3,e)=H_NN(3,7,e);
    H_NN(7,4,e)=H_NN(4,7,e);
    H_NN(7,5,e)=H_NN(5,7,e);    
    H_NN(7,6,e)=H_NN(6,7,e);
    H_NN(7,7,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,1)*J(e);
    H_NN(7,8,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,2)*J(e);
    H_NN(7,9,e)=1/L_mo*1/12*Nml_Vct(e,1)*Nml_Vct(e,3)*J(e);
    
    H_NN(8,1,e)=H_NN(1,8,e);
    H_NN(8,2,e)=H_NN(2,8,e);
    H_NN(8,3,e)=H_NN(3,8,e);
    H_NN(8,4,e)=H_NN(4,8,e);
    H_NN(8,5,e)=H_NN(5,8,e);    
    H_NN(8,6,e)=H_NN(6,8,e);    
    H_NN(8,7,e)=H_NN(7,8,e);  
    H_NN(8,8,e)=1/L_mo*1/12*Nml_Vct(e,2)*Nml_Vct(e,2)*J(e);
    H_NN(8,9,e)=1/L_mo*1/12*Nml_Vct(e,2)*Nml_Vct(e,3)*J(e);
    
    H_NN(9,1,e)=H_NN(1,9,e);
    H_NN(9,2,e)=H_NN(2,9,e);
    H_NN(9,3,e)=H_NN(3,9,e);
    H_NN(9,4,e)=H_NN(4,9,e);
    H_NN(9,5,e)=H_NN(5,9,e);    
    H_NN(9,6,e)=H_NN(6,9,e);    
    H_NN(9,7,e)=H_NN(7,9,e);
    H_NN(9,8,e)=H_NN(8,9,e);
    H_NN(9,9,e)=1/L_mo*1/12*Nml_Vct(e,3)*Nml_Vct(e,3)*J(e);
    
    for a=1:9
        for b=1:9
            % Assembly global H matrix by LM array
            P=LM(a,e); Q=LM(b,e);
            if (P>0) & (Q>0) 
                H(P,Q)=H_NN(a,b,e)+H(P,Q);
            end
        end
    end
end
% H_NN(:,:,1:20)

% prevent singular
% add small mass eps

eps=10^(-4);
 for i=1:B
     H(i,i)=H(i,i)+eps;
 end


disp('H OK!')

