function [vnodal] = fea(node,elem,mg,mp,M,gamma,ki)
%[vnodal] = fea(node,elem,m,M,gamma,ki)
%
%  vnodal: The nodal velocities
%    node: The nodal positions
%    elem: The element connectivities
%     m,M: Mobilities
%   gamma: Surface tension
%      ki: Used for energy anisotropy

%Sometimes this error will print. Use this command to turn it off
warning off MATLAB:nearlySingularMatrix;
%%[node,elem,mg,mp,M,gamma,ki,time,dt,solvertype,acc,mindt,maxdt,minel,maxel] = readdata('inputfile.dat');
%Determine the viscosity matrix
H=Hmatrix(node,elem,mg,mp,M);


%Determine the force vector
f=fvector(node,elem,gamma,ki);

%Solve the linear problem
q=H\f;


%Find the number of nodes
nnode=size(node,1);
%Find the number of elements
nel=size(elem,1);

%Extract the nodal velocities
vnodal=zeros(nnode,2);
% for i=1:nnode
%     if i==1
%         vnodal(i,1)=q(1);
%         vnodal(i,2)=q(2);
%     elseif i==2
%         vnodal(i,1)=q(4);
%         vnodal(i,2)=q(5);
%     else
%         vnodal(i,1)=q(4*(i-1));
%         vnodal(i,2)=q(4*(i-1)+1);
%     end
% end

for i=1:nel
    n1=elem(i,1);
    n2=elem(i,2);
    vnodal(n1,1)=q(3*n1-2);
    vnodal(n1,2)=q(3*n1-1);
    vnodal(n2,1)=q(3*n2-2);
    vnodal(n2,2)=q(3*n2-1);    
end