Uniaxial stretching solution not uniform in FEM code

I am trapped here for a long time. I wrote a toy Matlab FEM code. I want to run the follow simulation.

Mesh

Suppose we have a cube, and we divide it into subcube along $x,y,z$ axis, then each subcube be divide into 5 tetrahedral elements. The cube along $x,y,z$ axis is divide into $n_x,n_y,n_z$ node, respectively. So there are $n_x times n_y times n_z$ nodes in the grid. Each node in grid can be index with $node(i,j,k)$

Boundary conditions

• The cube’s $xoy$ plane has fixed the $z$ DOF, cube’s $xoz$ plane has fixed $y$ DOF, and cube’s $yoz$ plane has fixed $x$ DOF.

• The cube’s top plane has a uniform pressure applied (tensile).

So the cube is under uniaxial (tensile) stress.

Problems

I think under uniaxial stretching, the displacement along the z-axis of nodes having the same $z$ coordinate should be equal. That is $u_z(x, y, const) = const$.

However, my code doesn’t give me this solution. I do not know it is my code’s bug or it is the truth at they are not equal.

Main code

Actually, I already check the code with examples from some books. I already checked the element matrix and matrix assemble. My code can obtain the same solution with a standard book. Maybe there is something wrong with my loading. But I can not find a clue.

The full code and checking can be found in this repository

clc
clear
close all

num_x_node=4;
num_y_node=4;
num_z_node=4;

x_range=10;
y_range=10;
z_range=10;

grid= CubeDomainTetGrid(num_x_node, num_y_node, num_z_node, x_range, y_range, z_range);
node_coordinate_table=grid.nodeCoordinateTable();
element_node_table=grid.elementNodeTable();

dof_per_node=3;
num_node = grid.numNode();
num_dof = num_node * dof_per_node;

is_constrain = zeros(num_dof,1);

load_value = zeros(num_dof,1);
is_load = zeros(num_dof,1);

% assemble global stiffness matrix
element_stiffness_matrix_array = calculateElementStiffnessMatrixArray(node_coordinate_table, element_node_table);
K = assembleGlobalMatrix(node_coordinate_table, element_node_table, dof_per_node, element_stiffness_matrix_array);

%% constrain
% all node at x-y plane is fix z dof to zero
for i=1:1:num_x_node
    for j=1:1:num_y_node
        for k=1
            node_index = grid.nodeIndex(i,j,k);
            for dof_index = 3
                dof_global_index = (node_index - 1) * dof_per_node + dof_index;
                is_constrain(dof_global_index)=1;
            end
        end
    end
end
% all node at y-z plane is fix x dof to zero
for i=1:1:1
    for j=1:1:num_y_node
        for k=1:1:num_z_node
            node_index = grid.nodeIndex(i,j,k);
            for dof_index = 1
                dof_global_index = (node_index - 1) * dof_per_node + dof_index;
                is_constrain(dof_global_index)=1;
            end
        end
    end
end
% all node at x-z plane is fix y dof to zero
for i=1:1:num_x_node
    for j=1
        for k=1:1:num_z_node
            node_index = grid.nodeIndex(i,j,k);
            for dof_index = 2
                dof_global_index = (node_index - 1) * dof_per_node + dof_index;
                is_constrain(dof_global_index)=1;
            end
        end
    end
end
clear node_index;
clear dof_index;
clear dof_global_index;

%% load
% node(:,:,num_z_node) is apply pressure
pressure = 6;

% tranvers cell on surface and accumulate pressure
% the face see from top like
% the node force can be calulate by accumulate each small triangle
%  A y
%  |
%  |---> x
%  
%  -------------
%  |  |  |  |
%  |  |  |  |
%  |  |  |  |     
%  -------------
%  |  |  |  |
%  |  |  |  |
%  |  |  |  |
%  -------------
%
%      1   2
%      -----
%      |  |
%      |  |
%      |  |
%      -----
%      3   4

for i=1:1:num_x_node-1
    for j=1:1:num_y_node-1
        % calculate force on each sub cube's top surface
        dx= x_range / (num_x_node - 1);
        dy= y_range / (num_y_node - 1);
        area_triangle= dx*dy/2;
        force_on_triangle=area_triangle*pressure;
        num_node_per_triangle = 3;
        f_on_node = force_on_triangle/num_node_per_triangle;
        node_1=grid.nodeIndex(i,j+1,num_z_node);
        node_2=grid.nodeIndex(i+1,j+1,num_z_node);
        node_3=grid.nodeIndex(i,j,num_z_node);
        node_4=grid.nodeIndex(i+1,j,num_z_node);
        % lower triangle
        for node_index = [node_1 node_3 node_4]
            dof_index = 3;
            dof_global_index = (node_index - 1) * dof_per_node + dof_index;
            load_value(dof_global_index)=load_value(dof_global_index)+f_on_node;
        end
        % upper triangle
        for node_index = [node_1 node_2 node_4]
            dof_index = 3;
            dof_global_index = (node_index - 1) * dof_per_node + dof_index;
            load_value(dof_global_index)=load_value(dof_global_index)+f_on_node;
        end
        clear node_1 node_2 node_3 node_4
    end
end

%% apply constrain to global stiffness matrix
P=load_value;
for i=1:1:num_dof
    if is_constrain(i)
        K(i,:)=0;
        K(:,i)=0;
        K(i,i)=1;
        P(i)=0;
    end
end

%% solution
U=KP;

displacement=(reshape(U,dof_per_node,num_node))';
new_node_coordinate_table=node_coordinate_table+displacement*1e8;

figure
hold on
plotTetGrid(node_coordinate_table,element_node_table);
scatter3(new_node_coordinate_table(:,1),new_node_coordinate_table(:,2),new_node_coordinate_table(:,3));
axis equal

x_displacement=zeros(num_node,1);
y_displacement=zeros(num_node,1);
z_displacement=zeros(num_node,1);
for i=1:1:num_node
    x_displacement(i) = U((i - 1)*dof_per_node + 1);
    y_displacement(i) = U((i - 1)*dof_per_node + 2);
    z_displacement(i) = U((i - 1)*dof_per_node + 3);
end
figure
plot(z_displacement);
xlabel('node index');
ylabel('z displacement');
```

*node, nset=all_nodes
  1,                0,                0,                0
  2,                0,               10,                0
  3,               10,                0,                0
  4,               10,               10,                0
  5,                0,                0,               10
  6,                0,               10,               10
  7,               10,                0,               10
  8,               10,               10,               10
  9,                0,                0,               20
 10,                0,               10,               20
 11,               10,                0,               20
 12,               10,               10,               20
*element, type=c3d4, elset=all_elements
  1,   7,  10,  11,   8
  2,  12,  11,  10,   8
  3,   6,   3,   2,   4
  4,   6,   7,   3,   4
  5,   7,  10,   9,  11
  6,   6,  10,   7,   8
  7,   6,   8,   7,   4
  8,   2,   5,   1,   3
  9,   6,   5,   2,   3
 10,   6,   7,   5,   3
 11,   6,   9,   5,   7
 12,   6,  10,   9,   7
*material, name=the_mat
*elastic
         100.,  .3
*solid section, elset=all_elements, material=the_mat
1.0
*nset,nset=x0
1,2,5,6,9,10,
*nset,nset=y0
1,3,5,7,9,11,
*nset,nset=z0
1,2,3,4,
*nset,nset=top 
9,10,11,12,
*surface, TYPE=CUTTING SURFACE, DEFINITION=COORDINATES, name=top_surf
0.,0.,20.,   0.,0.,1.
all_elements
*Boundary
x0, 1, 1
y0, 2,2
z0, 3, 3
*Step, name=Step-1, nlgeom=NO
*Static
1., 1., 1e-05, 1.
*Dsload
**Surf-1, P, -1.
top_surf,P,-1.
*node print
u
*el print
S
*End Step

*Heading
** Job name: Job-test Model name: Job_my
** Generated by: Abaqus/CAE 2016
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=PART-1
*Node
      1,           0.,           0.,           0.
      2,           0.,          10.,           0.
      3,          10.,           0.,           0.
      4,          10.,          10.,           0.
      5,           0.,           0.,          10.
      6,           0.,          10.,          10.
      7,          10.,           0.,          10.
      8,          10.,          10.,          10.
      9,           0.,           0.,          20.
     10,           0.,          10.,          20.
     11,          10.,           0.,          20.
     12,          10.,          10.,          20.
*Element, type=C3D4
 1,  1,  3,  4,  7
 2,  1,  4,  2,  6
 3,  7,  5,  6,  1
 4,  7,  6,  8,  4
 5,  1,  7,  4,  6
 6,  5,  7,  8, 11
 7,  5,  8,  6, 10
 8, 11,  9, 10,  5
 9, 11, 10, 12,  8
10,  5, 11,  8, 10
*Elset, elset=Set-1, generate
  1,  10,   1
** Section: Section-1
*Solid Section, elset=Set-1, material=MATERIAL-1
,
*End Part
**  
**
** ASSEMBLY
**
*Assembly, name=Assembly
**  
*Instance, name=PART-1-1, part=PART-1
*End Instance
**  
*Nset, nset=Set-1, instance=PART-1-1, generate
 1,  4,  1
*Nset, nset=Set-2, instance=PART-1-1, generate
  1,  11,   2
*Nset, nset=Set-3, instance=PART-1-1
  1,  2,  5,  6,  9, 10
*Elset, elset=_Surf-1_S1, internal, instance=PART-1-1
 8, 9
*Surface, type=ELEMENT, name=Surf-1
_Surf-1_S1, S1
*End Assembly
** 
** MATERIALS
** 
*Material, name=MATERIAL-1
*Elastic
100.,0.
** 
** BOUNDARY CONDITIONS
** 
** Name: BC-1 Type: Displacement/Rotation
*Boundary
Set-1, 3, 3
** Name: BC-2 Type: Displacement/Rotation
*Boundary
Set-2, 2, 2
** Name: BC-3 Type: Displacement/Rotation
*Boundary
Set-3, 1, 1
** ----------------------------------------------------------------
** 
** STEP: Step-1
** 
*Step, name=Step-1, nlgeom=NO
*Static
1., 1., 1e-05, 1.
** 
** LOADS
** 
** Name: Load-1   Type: Pressure
*Dsload
Surf-1, P, -1.
** 
** OUTPUT REQUESTS
** 
*Restart, write, frequency=0
** 
** FIELD OUTPUT: F-Output-1
** 
*Output, field, variable=PRESELECT
** 
** HISTORY OUTPUT: H-Output-1
** 
*Output, history, variable=PRESELECT
*End Step