How to solve a singular matrix problem?

I am trying to write a code for 2d quadrature element using finite element method to solve the element stiffness matrix and force vector. However, my force vector has ended up with a singular matrix. How can I solve this problem? Below is the code I wrote.

%Mesh Information 
degP=1; %Number of nodes 
n=4; %Number of nodes per element
dof=1; %Number of dof per node 
ndof_elem=n*dof; %Number of dofs per element 
%Local node numbering scheme 
localNodsOnEdge = [1,2;2,3;3,4;4,1]; 
%Natural B.C.S for element 1 at each edge 
bcNv = {@(x,y)0,@(x,y)0,@(x,y)0,@(x,y)0}; %h = 0 at the edges of element 1 
%Material property: Thermal Conducivity (constant and isotropic) 
kij = @(x,y)(50+0.0*x+0.0*y);
%heat input
Q = @(x,y) sin(pi*x/2)*sin(pi*y/2);
%Gauss points and weights 
%Gauss Point in 2D
gp2d = [-1/sqrt(3.0)  -1/sqrt(3.0);
    1/sqrt(3.0)  -1/sqrt(3.0);
    1/sqrt(3.0)   1/sqrt(3.0);
    -1/sqrt(3.0)   1/sqrt(3.0)];
%Corresponding w is the weight for 2nd degree polynomial
w2d = [1.0;
    1.0;
    1.0;
    1.0];
%Shape functions and their derivatives 
%Evaluate shape functions at Gauss points for 2D case 
N2d = shapeFn2d(gp2d); %gp is the gauss points for 2nd degree polynomial  

%Evaluate derivative of the shape functions for 2D case
dN2d = dshapeFn(gp2d);
%Coordinate mesh information
conM = [1,2,3,4]; %Element 1 Connductivity matrix
corM = [0.5 0.5;1 0.5;1 1;0.5 1]; %corM = element 1 coordinate 
nElems = size(conM,1);
%Total nodes and degree of freedom, DOF
nNods = size(corM,1);
nDofs = nNods*dof;
%Global stiffness matrix and force vector 
K = zeros(nDofs);
F = zeros(nDofs,1);
for k=1:nElems
    ke=zeros(ndof_elem); %element stiffness matrix 
    fe=zeros(ndof_elem,1); %element force vector 
    cornerNods=localNodsOnEdge(:,1);
    %X and Y coordinates of the corner nodes (linear map)
    cordXk = corM(conM(k, cornerNods),1);  %X-coordinate of Corner nodes
    cordYk = corM(conM(k, cornerNods),2);  %Y-coordinate of Corner nodes
    cordXY = [cordXk cordYk]';
    %For each gauss point 
    for p = 1:size(gp2d,1)
        xGp = N2d(p,:)*cordXk;
        yGp = N2d(p,:)*cordYk;
        %calculting Jacobian at the Gauss point 
        J = cordXY*dN2d(:,:,p)';
        detJ = det(J);
        %Derivatives of shape functions in physical element
        dN_=(J\eye(2))'*dN2d(:,:,p);
        
        %calculating conductivity and heat infput at (xp,yp)
        kijGp = kij(xGp, yGp);
        QGp = Q(xGp,yGp);
        
        %calculating element stiffness and force entries 
        for i = 1:ndof_elem  %dof per element 
            for j = 1:ndof_elem
                ke(i,j) = ke(i,j)+dN_(:,i)'*kijGp*dN_(:,j)*detJ*w2d(p);
            end
            fe(i) = fe(i)+QGp*N2d(p, i)*detJ*w2d(p);
        end
    end
    %Assembly
    for i=1:ndof_elem
        ig=conM(k,i);
        for j=1:ndof_elem
            jg=conM(k,j);
            K(ig,jg)=K(ig,jg)+ke(i,j);
        end
        F(ig)=F(ig)+fe(i);
    end
end
K
F
%obtian the finite element solutions 
u=K\F

function dN=dshapeFn(gp)
%All shape functions are arranged row-wise
n_gp=size(gp,1);
dN=zeros(2,4,n_gp);
for k=1:n_gp
    %ksi=gp(k,1); eta=gp(k,2)
    dN(:,:,k)=[-(1/4)*(1-gp(k,2)) (1/4)*(1-gp(k,2))  (1/4)*(1+gp(k,2)) -(1/4)*(1+gp(k,2));
        -(1/4)*(1-gp(k,1)) -(1/4)*(1+gp(k,1)) (1/4)*(1+gp(k,1)) (1/4)*(1-gp(k,1))];
end
end 

function N=shapeFn2d(gp)
%ksi and eta are column vectors 
%ksi=gp(:,1);
%eta=gp(:,2);
N=[(1/4)*(1-gp(:,1)).*(1-gp(:,2)) (1/4)*(1+gp(:,1)).*(1-gp(:,2))...        
    (1/4)*(1+gp(:,1)).*(1+gp(:,2)) (1/4)*(1-gp(:,1)).*(1+gp(:,2))];
end

Answer 1

[John D'Errico]

Probable user error.

How can you solve it? Almost always with a FEM problem, a singularity arises when you have made a mistake. This is usually a great clue. :)

What kind of mistake are we talking about? There are several common ones I can think of.

First, you may have failed to fully specify sufficient boundary conditions, so that the solution is ambiguous. Essentially there are infinitely many solutions. What do I mean here? I'll offer an example. Suppose I build a model of a bridge. So a simple truss problem. But I need to specify that the x and y displacement of at least some point in the truss is fixed. If not, then you can freely translate the truss to any location, without paying any energy penalty. What will you see then? A SINGULAR MATRIX! Essentially, this tells you that there are infinitely many solutions, and all are equally valid in terms of energy penalty.

There are other ways a FEM can fail due to a singularity. For example, if you allow it to freely rotate around a point, with no energy penalty, then again, you will see a singular matrix.

What else can you do wrong? You might have duplicate nodes. So two nodes at the same location. That will cause zero volume elements, or zero length elements, in a truss model, for example. Again, you would expect a singularity to arise.

As you can see, these most common errors one can make are really an error in constructing the model. What can you do about them? DON'T MAKE AN ERROR. Really, what more can I say? I would look very carefully at what you have specified. Did you make a mistake?

[Walter Roberson]

From what I have observed, it is common to forget to prevent z rotations -- to lock in x and y on both sides but forget that the span could spin.

(And remember that you do not want only one point being the barrier to rotation: that puts a lot of strain on that one point.)

[John D'Errico]

Yes. Back when I was an engineer (not the kind who make the trains run either) it was rotation that was easy to miss. But anything that allows the system to move without deformation is a cause for singularity, because then there are infinitely many solutions.