# Why the solid FEM problem can not be solved after constraining 3 degrees of freedom?

I write a simple MATLAB code for solving solid FEM problem.

The problem looks like that

(1)      (2)
x-------x
|     / |
|   /   |
| /     |
x-------x
(3)     (4)


Each node has two degrees of freedom(DOF), which stand for $$x$$ and $$y$$ displacements.

# Boundary condition

I want only constraint node 1's $$x,y$$ DOF and node 2's \$x DOF, that is constraining DOFs 1,2,3 to zero.

All nodes have a node force equal to 1.

# Result

I found that the global stiffness matrix has a rank 7 but not the rank 8, which leads to the problem can not be solved.

# Problem

In the 2D solid FEM case, the 3 DOFs are needed to be constrained for limiting rigid body translation along $$x,y$$ axis, and the rotaion. But after I constrain 3 DOFs, the problem still can not be solved. Why?

# Appendix

The matlab code looks like

clc
clear
close all
%% model
x_range=10;
y_range=10;
num_x_node=2;
num_y_node=2;
num_node=num_x_node*num_y_node;
DOF=2;
%% generate grid
node_coordinate_table=calculateNodeCoordinateTableForRectangleDomain(num_x_node,num_y_node,x_range,y_range);
element_node_table=calculateElementNodeTableForRectangleDomainTriangleGrid(num_x_node,num_y_node);
%% calculate element matrix
K=calcGlobalMatrix(node_coordinate_table,element_node_table,DOF);
P=ones(num_node*DOF,1);
%% constrain
is_fix_dof=zeros(num_node*DOF,1);
is_fix_dof(1)=1;
is_fix_dof(2)=1;
is_fix_dof(3)=1;
%% modify discrete equation by constrain (cross zero center one method)
for i=1:1:num_node*DOF
if is_fix_dof(i)==1
K(i,:)=0;
K(:,i)=0;
K(i,i)=1;
P(i)=0;
end
end
%% solution
x=K\P
%% plot grid
figure
set(gcf,'position',[50,50,1000,800])
hold on
plotNode(node_coordinate_table);
plotElement(element_node_table,node_coordinate_table);

%% function
function K=calcGlobalMatrix(node_coordinate_table,element_node_table,DOF)
num_node=size(node_coordinate_table,1);
K=zeros(DOF*num_node);
element_matrix_array=calculateElementMatrixArray(node_coordinate_table,element_node_table);
num_element=size(element_node_table,1);
for i=1:1:num_element
K0=element_matrix_array{i};
node_global_index_array=element_node_table(i,:);
for node_i=1:1:3
for node_j=1:1:3
for i1=1:1:DOF
for j1=1:1:DOF
node_global_index_i=node_global_index_array(node_i);
node_global_index_j=node_global_index_array(node_j);
row=(node_i-1)*DOF+i1;
col=(node_j-1)*DOF+j1;
ROW=(node_global_index_i-1)*DOF+i1;
COL=(node_global_index_j-1)*DOF+j1;
K(ROW,COL)=K(ROW,COL)+K0(row,col);
end
end
end
end
end

end

function [element_matrix_array]=calculateElementMatrixArray(node_coordinate_table,element_node_table)
num_element=size(element_node_table,1);
element_matrix_array=cell(num_element,1);
for i=1:1:num_element
node_a_id=element_node_table(i,1);
node_b_id=element_node_table(i,2);
node_c_id=element_node_table(i,3);
x=zeros(3,1);
x(1)=node_coordinate_table(node_a_id,1);
x(2)=node_coordinate_table(node_b_id,1);
x(3)=node_coordinate_table(node_c_id,1);
y=zeros(3,1);
y(1)=node_coordinate_table(node_a_id,2);
y(2)=node_coordinate_table(node_b_id,2);
y(3)=node_coordinate_table(node_c_id,2);
element_matrix_array{i}=calculateElementMatrixByNodeCoordinate(x,y);
end

end

function plotNode(node_coordinate_table)
x=node_coordinate_table(:,1);
y=node_coordinate_table(:,2);
plot(x,y,'o');
num_node=size(node_coordinate_table,1);
for i=1:1:num_node
node_x=x(i);
node_y=y(i);
text(node_x,node_y,num2str(i));
end

end

function plotElement(element_node_table,node_coordinate)
num_element=size(element_node_table,1);
for i=1:1:num_element
node_a_id=element_node_table(i,1);
node_b_id=element_node_table(i,2);
node_c_id=element_node_table(i,3);
x_a=node_coordinate(node_a_id,1);
y_a=node_coordinate(node_a_id,2);
x_b=node_coordinate(node_b_id,1);
y_b=node_coordinate(node_b_id,2);
x_c=node_coordinate(node_c_id,1);
y_c=node_coordinate(node_c_id,2);
plot([x_a,x_b,x_c,x_a],[y_a,y_b,y_c,y_a],'k-');
x_o=(x_a+x_b+x_c)/3;
y_o=(y_a+y_b+y_c)/3;
text(x_o,y_o,[num2str(i)]);% -1 for plot element index
end
end

function node_coordinate_table=calculateNodeCoordinateTableForRectangleDomain(num_x_node,num_y_node,x_range,y_range)
num_node=num_x_node*num_y_node;
num_x_element=num_x_node-1;
num_y_element=num_y_node-1;
node_coordinate_table=zeros(num_node,2);
n=1;
dx=x_range/num_x_element;
dy=y_range/num_y_element;
for j=1:1:num_y_node
for i=1:1:num_x_node
x = (i-1) * dx;
y = y_range - (j-1) * dy;
node_coordinate_table(n,1)=x;
node_coordinate_table(n,2)=y;
n=n+1;
end
end
end

function element_node_table=calculateElementNodeTableForRectangleDomainTriangleGrid(num_x_node,num_y_node)
% element node table generate is order in matlab format, that is, id start from 1
num_x_element=num_x_node-1;
num_y_element=num_y_node-1;
num_element=num_x_element*num_y_element*2;
element_node_table = zeros(num_element,3);
n=1;
for j=1:1:num_y_element
for i=1:1:num_x_element
%  a--c
%  | /|
%  |/ |
%  b--d
node_a_id=i+(j-1)*num_x_node;
node_b_id=node_a_id+num_x_node;
node_c_id=node_a_id+1;
node_d_id=node_b_id+1;
% upper element abc
element_node_table(n,1)=node_a_id;
element_node_table(n,2)=node_b_id;
element_node_table(n,3)=node_c_id;
n=n+1;
% lower element cbd
element_node_table(n,1)=node_c_id;
element_node_table(n,2)=node_b_id;
element_node_table(n,3)=node_d_id;
n=n+1;
end
end
end

function K0=calculateElementMatrixByNodeCoordinate(x,y)
% reference
% finite element method, wang xucheng chapter 2

% model
E_0=100;
nu_0=0;
Thickness=1;
% calculate element matrix
a=zeros(3,1);
b=zeros(3,1);
c=zeros(3,1);
for i=1:1:3
if(i==1);j=2;m=3;end
if(i==2);j=3;m=1;end
if(i==3);j=1;m=2;end
a(i)=x(j)*y(m)-x(m)*y(j);
b(i)=y(j)-y(m);
c(i)=-x(j)+x(m);
end
A=1/2*(b(i)*c(j)-b(j)*c(i));
K0=zeros(6,6);
for r=1:1:3
for s=1:1:3
K0(2*r-1,2*s-1)=b(r)*b(s)+(1-nu_0)/2*c(r)*c(s);
K0(2*r-1,2*s)=nu_0*b(r)*c(s)+(1-nu_0)/2*c(r)*b(s);
K0(2*r,2*s-1)=nu_0*c(r)*b(s)+(1-nu_0)/2*b(r)*c(s);
K0(2*r,2*s)=c(r)*c(s)+(1-nu_0)/2*b(r)*b(s);
end
end
K0=K0*E_0*Thickness/(4*A*(1-nu_0^2));
end
$$$$
`