Single hexahedral element stiffness matrix problem, help me find the mistake

Attached below is some code I wrote to solve a basic problem: finding the node displacements of a cube with two vertices constrained (vertices 6 and 7 with coordinates [1,1,1], [1,-1,1]), and the rest of the vertices loaded with point loads in the y direction. What I expect is for the displacement at vertex 0, with coordinates [-1,-1,-1], to be strictly inside the plane spanned by [0,1,0], [1,0,1] because of the symmetry of the problem. What I get is a displacement that is not on that plane, with components [-0.00289863, 0.07050896, -0.02001042].

The way I solve the problem is by calculating the stiffness matrix with Gaussian integration. Then I remove the rows and columns in the matrix that correspond to the constrained vertices (actually I set them to 0, and the element on the diagonal to 1, but it should give the same results, right?). And then I solve the equation kd = f where k is the stiffness matrix, d is the displacement vector (24 dimensions) and f is the nodal forces vector.

I tested the shape function derivatives (test code included) by numerically calculating the derivatives and comparing with the ones returned by the function. I also tested the strain matrix by calculating the strain numerically. These are not to blame.

Please, help me find my mistake.

Edit: The problem is not well-formed, because there is nothing preventing the cube from rotating around the vertices 6-7 axis. However, why is the determinant of the stiffness matrix not 0? Shouldn't the matrix be singular?

import numpy as np

from math import sqrt

_samplePoints = [[0.0],
[-sqrt(1./3),sqrt(1./3)],
[-sqrt(0.6),0.0,sqrt(0.6)],
[-0.861136, -0.339981, 0.339981, 0.861136],
[-0.906180, -0.538469, 0, 0.538469, 0.906180],
[-0.932470, -0.661209, -0.238619, 0.238619, 0.661209, 0.932470]]
_coeffs = [[2.0],
[1.0,1.0],
[0.555556,0.888889,0.555556],
[0.347855, 0.652145, 0.652145, 0.347855],
[0.236927, 0.478629, 0.568889,0.478629, 0.236927],
[0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324]]

return _samplePoints[numPoints-1],_coeffs[numPoints-1]

_naturalCoordinatesVertices = np.array([[-1,-1,-1],[-1,1,-1],[-1,1,1],[-1,-1,1],[1,-1,-1],[1,1,-1],[1,1,1],[1,-1,1]])

def ShapeFuncVec(ksiEtaZeta):
return 0.125 * np.prod(1 + ksiEtaZeta.reshape(1,3) * _naturalCoordinatesVertices,axis=1)

def DShapeFuncVec_DNaturalCoordinates(ksiEtaZeta):
return 0.125 * _naturalCoordinatesVertices[:,[0,1,2]]*(1+_naturalCoordinatesVertices[:,[1,2,0]]*ksiEtaZeta[[1,2,0]])*(1+_naturalCoordinatesVertices[:,[2,0,1]]*ksiEtaZeta[[2,0,1]])

def CalcJacobian(ksiEtaZeta,vertices):
return DShapeFuncVec_DNaturalCoordinates(ksiEtaZeta).transpose().dot(vertices)

def DShapeFuncVec_DXYZ(ksiEtaZeta,vertices):
dShapeFuncVec_DKsiEtaZeta = DShapeFuncVec_DNaturalCoordinates(ksiEtaZeta)
J = dShapeFuncVec_DKsiEtaZeta.transpose().dot(vertices)
Jinv = np.linalg.inv(J)
dShapeFuncVec_dXYZ = dShapeFuncVec_DKsiEtaZeta.dot(Jinv.transpose())
return dShapeFuncVec_dXYZ, J

def CalcStrainMatrix(ksiEtaZeta, vertices):
dN_dxyz,J = DShapeFuncVec_DXYZ(ksiEtaZeta,vertices)
assert(np.linalg.det(J)>0)
B = np.zeros((6,8*3))
for i in range(8):
Bi = np.zeros((6,3))
dNi_dxyz = dN_dxyz[i,:]
Bi[[0,1,2],[0,1,2]] = dNi_dxyz
Bi[[3,3,4,4,5,5],[1,2,0,2,0,1]] = dNi_dxyz[[2,1,2,0,1,0]]
B[:,(i*3):(i*3+3)] = Bi
return B,J

def GetMaterialTensor():
youngsModulus = 50.0
poissonRatio = 0.4
denom = ((1-2*poissonRatio)*(1+poissonRatio))
c11 = youngsModulus * (1-poissonRatio) / denom
c12 = youngsModulus * poissonRatio / denom
materialTensor = np.zeros((6,6))
materialTensor[:3,:3] = c12
materialTensor[[0,1,2],[0,1,2]] = c11
materialTensor[[3,4,5],[3,4,5]] = (c11-c12)/2
return materialTensor

materialTensor = GetMaterialTensor()

import itertools

def CalcStiffnessMatrix(vertices, materialTensor):
n=3
#perform gauss integration with n points
Ke = np.zeros((24,24),dtype=float)
count=0
for i,j,k in itertools.product(range(n),repeat=3):
count+=1
ksiEtaZeta = np.array([samplePoints[i],samplePoints[j],samplePoints[k]])
coeff = coeffs[i]*coeffs[j]*coeffs[k]
B,J = CalcStrainMatrix(ksiEtaZeta,vertices)
Ke += np.linalg.det(J)*(B.transpose().dot(materialTensor).dot(B)) * coeff
assert(count == n**3)
return Ke

def TestDerivatives():
ksiEtaZeta_sample = np.array([0.4,0.3,0.2])

print(ShapeFuncVec(ksiEtaZeta_sample))
print(sum(ShapeFuncVec(ksiEtaZeta_sample)))

print(DShapeFuncVec_DNaturalCoordinates(ksiEtaZeta_sample))

print(CalcJacobian(ksiEtaZeta_sample,vertices))

dksiEtaZeta = np.array([0.000001,-0.00002,0.000003])

xyz_sample = ShapeFuncVec(ksiEtaZeta_sample).dot(vertices)
xyz_sample2 = ShapeFuncVec(ksiEtaZeta_sample + dksiEtaZeta).dot(vertices)
print(xyz_sample,xyz_sample2)

dxyz_sample = xyz_sample2 - xyz_sample
dxyz_sample_estimated = CalcJacobian(ksiEtaZeta_sample,vertices).transpose().dot(dksiEtaZeta)
print(dxyz_sample)
print(dxyz_sample_estimated)

assert(np.max(np.abs(dxyz_sample-dxyz_sample_estimated))/np.max(np.abs(dxyz_sample)) < 1e-5)

x1 = ShapeFuncVec(ksiEtaZeta_sample + dksiEtaZeta) - ShapeFuncVec(ksiEtaZeta_sample)
x2 = DShapeFuncVec_DXYZ(ksiEtaZeta_sample,vertices)[0].dot(dxyz_sample)
print(x1)
print(x2)

assert(np.max(np.abs(x1-x2))/np.max(np.abs(x1)) < 1e-5)

def TestStrain():
ksiEtaZeta = np.random.rand(3) *2-1
B,J = CalcStrainMatrix(ksiEtaZeta,vertices)
# print(J)

#displacements = np.array([[0.1,0.5,-0.1]] + [[0,0,0]]*7)
displacements = 0.2 * np.random.rand(8,3)-0.1
displacements_P = displacements.transpose().dot(ShapeFuncVec(ksiEtaZeta))

strain_fromB = B.dot(displacements.flatten())

strain_fromNumDiff = np.zeros((3,3))
dKezMag = 0.01
dKez = np.eye(3) * dKezMag
kezPdkez = np.vstack([ksiEtaZeta]*3) + dKez
for i in range(3):
thisDisplacement = displacements.transpose().dot(ShapeFuncVec(kezPdkez[i,:]))
dDisplacement = thisDisplacement - displacements_P
strain_fromNumDiff[:,i] = dDisplacement / dKezMag

strain_fromNumDiff_vec = np.array([strain_fromNumDiff[0,0],
strain_fromNumDiff[1,1],
strain_fromNumDiff[2,2],
strain_fromNumDiff[1,2]+strain_fromNumDiff[2,1],
strain_fromNumDiff[0,2]+strain_fromNumDiff[2,0],
strain_fromNumDiff[1,0]+strain_fromNumDiff[0,1]])

print(strain_fromB)
print(strain_fromNumDiff)
print(strain_fromNumDiff_vec)

assert(np.max(np.abs(strain_fromB-strain_fromNumDiff_vec)) < 1e-6)
print(np.max(np.abs(strain_fromB-strain_fromNumDiff_vec)))

vertices = _naturalCoordinatesVertices.copy()
Ke = CalcStiffnessMatrix(vertices,materialTensor)
constrainedVertices = np.vstack([[6,7]]*3)
constrainedVertices *= 3
constrainedVertices += np.arange(3).reshape(3,-1)
Ke[constrainedVertices,:]=0
Ke[:,constrainedVertices]=0
Ke[constrainedVertices,constrainedVertices]=1
print("det(Ke) = ",np.linalg.det(Ke))
print("eig(Ke) = ",np.linalg.eig(Ke)[0])
forces = np.zeros((8,3),dtype=float)
forces[:6,1]=0.1
forces=forces.flatten()
displacements = np.linalg.solve(Ke,forces).reshape(-1,3)
print()
print("result nodal displacements = ")
print(displacements)

print()
print("Tests:")

TestDerivatives()

TestStrain()


Output:

det(Ke) =  12856964.398889132
eig(Ke) =  [1.98207971e+02 5.97329352e+01 5.36174927e+01 4.56626826e+01
3.28979711e+01 3.02076122e+01 2.88176063e+01 2.70455110e+01
2.10812642e+01 1.73934001e+01 1.69113782e+01 1.29241062e+01
9.90047362e+00 9.82465058e+00 4.45269904e+00 9.04339639e-15
1.37800889e+00 1.37375076e+00 1.00000000e+00 1.00000000e+00
1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00]

result nodal displacements =
[[-0.00289863  0.07050896 -0.02001042]
[ 0.02001042  0.07050896  0.00289863]
[ 0.01565451  0.04365447  0.01244406]
[-0.01565451  0.04365447 -0.02955586]
[-0.01244406  0.04365447 -0.01565451]
[ 0.02955586  0.04365447  0.01565451]
[ 0.          0.          0.        ]
[ 0.          0.          0.        ]]

Tests:
[0.042 0.078 0.117 0.063 0.098 0.182 0.273 0.147]
1.0
[[-0.07   -0.06   -0.0525]
[-0.13    0.06   -0.0975]
[-0.195   0.09    0.0975]
[-0.105  -0.09    0.0525]
[ 0.07   -0.14   -0.1225]
[ 0.13    0.14   -0.2275]
[ 0.195   0.21    0.2275]
[ 0.105  -0.21    0.1225]]
[[ 1.00000000e+00 -1.38777878e-17  1.38777878e-17]
[ 0.00000000e+00  1.00000000e+00  0.00000000e+00]
[-1.38777878e-17  0.00000000e+00  1.00000000e+00]]
[0.4 0.3 0.2] [0.400001 0.29998  0.200003]
[ 1.e-06 -2.e-05  3.e-06]
[ 1.e-06 -2.e-05  3.e-06]
[ 9.72493763e-07 -1.62249301e-06 -1.70250199e-06  1.85250124e-06
2.50249124e-06 -3.35249199e-06 -3.32251301e-06  4.67251376e-06]
[ 9.7250e-07 -1.6225e-06 -1.7025e-06  1.8525e-06  2.5025e-06 -3.3525e-06
-3.3225e-06  4.6725e-06]
[-0.01654885 -0.00454239  0.00107933  0.01714707 -0.00310323 -0.09393567]
[[-0.01654885 -0.03663831 -0.00499938]
[-0.05729736 -0.00454239  0.00809046]
[ 0.00189615  0.00905662  0.00107933]]
[-0.01654885 -0.00454239  0.00107933  0.01714707 -0.00310323 -0.09393567]
7.927686285214008e-16

• Your stiffness matrix calculation is correct. I think the problem is you are not constraining the DOFs you think you are. Take a look at Ke after the lines that apply the constraints. Aug 19 '21 at 19:47

Your intuition is right, the stiffness matrix should be singular. Actually, as pointed by @BillGreene, it has an eigenvalue that is zero.

Ordering the nodes as follows (-1, -1, -1,), (1, -1, -1), (1, 1, 1), (-1, 1, -1), (-1, -1, 1,), (1, -1, 1), (1, 1, 1), (-1, 1, 1) and constraining the nodes 0 and 1, we get the following stiffness matrix (for $$\nu=4/10$$)

$$K/E = \begin{bmatrix} I_6 &0\\ 0 &K_2 \end{bmatrix}\, ,$$

with

$$K_2 = \left[\begin{array}{cccccccccccccccccc}\frac{40}{63} & \frac{25}{84} & - \frac{25}{84} & - \frac{25}{63} & \frac{5}{28} & - \frac{5}{28} & - \frac{10}{63} & - \frac{25}{168} & \frac{25}{168} & \frac{5}{126} & - \frac{5}{56} & \frac{5}{56} & \frac{25}{126} & \frac{25}{168} & \frac{5}{28} & - \frac{65}{252} & \frac{5}{56} & \frac{25}{84}\\\frac{25}{84} & \frac{40}{63} & - \frac{25}{84} & - \frac{5}{28} & \frac{25}{126} & - \frac{25}{168} & - \frac{25}{168} & - \frac{10}{63} & \frac{25}{168} & \frac{5}{56} & - \frac{65}{252} & \frac{25}{84} & \frac{25}{168} & \frac{25}{126} & \frac{5}{28} & - \frac{5}{56} & \frac{5}{126} & \frac{5}{56}\\- \frac{25}{84} & - \frac{25}{84} & \frac{40}{63} & \frac{5}{28} & - \frac{25}{168} & \frac{25}{126} & \frac{25}{168} & \frac{25}{168} & - \frac{10}{63} & - \frac{5}{56} & \frac{25}{84} & - \frac{65}{252} & - \frac{5}{28} & - \frac{5}{28} & - \frac{25}{63} & \frac{25}{84} & - \frac{5}{56} & - \frac{65}{252}\\- \frac{25}{63} & - \frac{5}{28} & \frac{5}{28} & \frac{40}{63} & - \frac{25}{84} & \frac{25}{84} & \frac{5}{126} & \frac{5}{56} & - \frac{5}{56} & - \frac{10}{63} & \frac{25}{168} & - \frac{25}{168} & - \frac{65}{252} & - \frac{5}{56} & - \frac{25}{84} & \frac{25}{126} & - \frac{25}{168} & - \frac{5}{28}\\\frac{5}{28} & \frac{25}{126} & - \frac{25}{168} & - \frac{25}{84} & \frac{40}{63} & - \frac{25}{84} & - \frac{5}{56} & - \frac{65}{252} & \frac{25}{84} & \frac{25}{168} & - \frac{10}{63} & \frac{25}{168} & \frac{5}{56} & \frac{5}{126} & \frac{5}{56} & - \frac{25}{168} & \frac{25}{126} & \frac{5}{28}\\- \frac{5}{28} & - \frac{25}{168} & \frac{25}{126} & \frac{25}{84} & - \frac{25}{84} & \frac{40}{63} & \frac{5}{56} & \frac{25}{84} & - \frac{65}{252} & - \frac{25}{168} & \frac{25}{168} & - \frac{10}{63} & - \frac{25}{84} & - \frac{5}{56} & - \frac{65}{252} & \frac{5}{28} & - \frac{5}{28} & - \frac{25}{63}\\- \frac{10}{63} & - \frac{25}{168} & \frac{25}{168} & \frac{5}{126} & - \frac{5}{56} & \frac{5}{56} & \frac{40}{63} & \frac{25}{84} & - \frac{25}{84} & - \frac{25}{63} & \frac{5}{28} & - \frac{5}{28} & - \frac{65}{252} & - \frac{25}{84} & - \frac{5}{56} & \frac{25}{126} & - \frac{5}{28} & - \frac{25}{168}\\- \frac{25}{168} & - \frac{10}{63} & \frac{25}{168} & \frac{5}{56} & - \frac{65}{252} & \frac{25}{84} & \frac{25}{84} & \frac{40}{63} & - \frac{25}{84} & - \frac{5}{28} & \frac{25}{126} & - \frac{25}{168} & - \frac{25}{84} & - \frac{65}{252} & - \frac{5}{56} & \frac{5}{28} & - \frac{25}{63} & - \frac{5}{28}\\\frac{25}{168} & \frac{25}{168} & - \frac{10}{63} & - \frac{5}{56} & \frac{25}{84} & - \frac{65}{252} & - \frac{25}{84} & - \frac{25}{84} & \frac{40}{63} & \frac{5}{28} & - \frac{25}{168} & \frac{25}{126} & \frac{5}{56} & \frac{5}{56} & \frac{5}{126} & - \frac{25}{168} & \frac{5}{28} & \frac{25}{126}\\\frac{5}{126} & \frac{5}{56} & - \frac{5}{56} & - \frac{10}{63} & \frac{25}{168} & - \frac{25}{168} & - \frac{25}{63} & - \frac{5}{28} & \frac{5}{28} & \frac{40}{63} & - \frac{25}{84} & \frac{25}{84} & \frac{25}{126} & \frac{5}{28} & \frac{25}{168} & - \frac{65}{252} & \frac{25}{84} & \frac{5}{56}\\- \frac{5}{56} & - \frac{65}{252} & \frac{25}{84} & \frac{25}{168} & - \frac{10}{63} & \frac{25}{168} & \frac{5}{28} & \frac{25}{126} & - \frac{25}{168} & - \frac{25}{84} & \frac{40}{63} & - \frac{25}{84} & - \frac{5}{28} & - \frac{25}{63} & - \frac{5}{28} & \frac{25}{84} & - \frac{65}{252} & - \frac{5}{56}\\\frac{5}{56} & \frac{25}{84} & - \frac{65}{252} & - \frac{25}{168} & \frac{25}{168} & - \frac{10}{63} & - \frac{5}{28} & - \frac{25}{168} & \frac{25}{126} & \frac{25}{84} & - \frac{25}{84} & \frac{40}{63} & \frac{25}{168} & \frac{5}{28} & \frac{25}{126} & - \frac{5}{56} & \frac{5}{56} & \frac{5}{126}\\\frac{25}{126} & \frac{25}{168} & - \frac{5}{28} & - \frac{65}{252} & \frac{5}{56} & - \frac{25}{84} & - \frac{65}{252} & - \frac{25}{84} & \frac{5}{56} & \frac{25}{126} & - \frac{5}{28} & \frac{25}{168} & \frac{40}{63} & \frac{25}{84} & \frac{25}{84} & - \frac{25}{63} & \frac{5}{28} & \frac{5}{28}\\\frac{25}{168} & \frac{25}{126} & - \frac{5}{28} & - \frac{5}{56} & \frac{5}{126} & - \frac{5}{56} & - \frac{25}{84} & - \frac{65}{252} & \frac{5}{56} & \frac{5}{28} & - \frac{25}{63} & \frac{5}{28} & \frac{25}{84} & \frac{40}{63} & \frac{25}{84} & - \frac{5}{28} & \frac{25}{126} & \frac{25}{168}\\\frac{5}{28} & \frac{5}{28} & - \frac{25}{63} & - \frac{25}{84} & \frac{5}{56} & - \frac{65}{252} & - \frac{5}{56} & - \frac{5}{56} & \frac{5}{126} & \frac{25}{168} & - \frac{5}{28} & \frac{25}{126} & \frac{25}{84} & \frac{25}{84} & \frac{40}{63} & - \frac{5}{28} & \frac{25}{168} & \frac{25}{126}\\- \frac{65}{252} & - \frac{5}{56} & \frac{25}{84} & \frac{25}{126} & - \frac{25}{168} & \frac{5}{28} & \frac{25}{126} & \frac{5}{28} & - \frac{25}{168} & - \frac{65}{252} & \frac{25}{84} & - \frac{5}{56} & - \frac{25}{63} & - \frac{5}{28} & - \frac{5}{28} & \frac{40}{63} & - \frac{25}{84} & - \frac{25}{84}\\\frac{5}{56} & \frac{5}{126} & - \frac{5}{56} & - \frac{25}{168} & \frac{25}{126} & - \frac{5}{28} & - \frac{5}{28} & - \frac{25}{63} & \frac{5}{28} & \frac{25}{84} & - \frac{65}{252} & \frac{5}{56} & \frac{5}{28} & \frac{25}{126} & \frac{25}{168} & - \frac{25}{84} & \frac{40}{63} & \frac{25}{84}\\\frac{25}{84} & \frac{5}{56} & - \frac{65}{252} & - \frac{5}{28} & \frac{5}{28} & - \frac{25}{63} & - \frac{25}{168} & - \frac{5}{28} & \frac{25}{126} & \frac{5}{56} & - \frac{5}{56} & \frac{5}{126} & \frac{5}{28} & \frac{25}{168} & \frac{25}{126} & - \frac{25}{84} & \frac{25}{84} & \frac{40}{63}\end{array}\right]\, .$$

For a load vector aligned with the edge, i.e,

$$f^T = \left[\begin{array}{cccccccccccccccccccccccc}0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0\end{array}\right]\, ,$$

the solution is the following

$$u = \left( 0, \ 0, \ 0, \ 0, \ 0, \ 0, \ \frac{2401}{110}, \ - \frac{861}{110}, \ \tau_{0} - \frac{357}{22}, \ \frac{2401}{110}, \ \frac{861}{110}, \ \tau_{0} + \frac{105}{22}, \ \frac{2401}{110}, \ \frac{357}{22} - \tau_{0}, \ \frac{861}{110}, \ \frac{2401}{110}, \ - \tau_{0} - \frac{105}{22}, \ - \frac{861}{110}, \ \frac{1939}{55}, \ - \tau_{0}, \ \tau_{0} - \frac{126}{11}, \ \frac{1939}{55}, \ \frac{126}{11} - \tau_{0}, \ \tau_{0}\right)\, ,$$

with $$\tau_0$$ a free parameter. Notice that

$$\tau_{0}\left( 0, \ 0, \ 0, \ 0, \ 0, \ 0, \ 0, \ 0, \ 1, \ 0, \ 0, \ 1, \ 0, \ - 1, \ 0, \ 0, \ -1, \ 0, \ 0, \ - 1, \ 1, \ 0, \ - 1, \ 1\right)\, ,$$

is the eigenspace associated with the null eigenvalue.

• Very good answer! Aug 19 '21 at 22:23
• I'll add to that, that there was a tiny eigen value hiding in the print-out: 9.04339639e-15. Because the np.linalg.eig didn't actually order the eigenvalues as I expected, I didn't notice it. So I guess that the matrix was singular in my solution as well. What's funny is that np.linalg.det actually returned a fairly large number, so I think there is some numerical instability in numpy that leads to the determinant being large while there is an almost zero eigen value. Aug 20 '21 at 0:10

Each node has only 3 degrees of freedom: UX, UY and UZ. As you constrained nodes 6 and 7 in all three directions, there is no rigid body motion. Your problem is well-posed. This is why the stiffness matrix should not be singular. In other words, its determinant should be non-zero.

Since there are no rotational DOFs, you should not observe any rotation, unless your loading conditions force such a scenario. But, your load is only in Y-direction; so, there should not be any rotation. You should observe bending in Y-direction.

As I am not an expert user of Python, your code is too cryptic for me to follow. Nevertheless, if your stiffness is not symmetric, then there is a bug in your code. Likely to be in the portion of the code where you calculate the derivatives of the shape functions. Why do you use abs when you add the contribution to Ke? You should not do that.

Instead of computing the stiffness matrix first, calculate the mass matrix and volume first. This will help in finding the bug.

• 1. Does the rotation DOF need to be specified explicitly? If it's a free body, wouldn't it rotate if given torque around the y axis? 2. abs was used because I was using the jacobian to integrate volume and volume can't be negative. But I removed the abs in any case. Also, there's an assertion elsewhere that det(J)>0. 3. I numerically tested the derivatives of the shape functions and the strain matrix. They correctly predict the strain as a result of a change of the vertices coordinates. Aug 18 '21 at 13:36
• I just added the code to test the shape function derivatives and the strain matrix calculation. Aug 18 '21 at 13:50
• Continuum elements do not have rotation DOFs. So, no need to specify them. I think it would be helpful if you can post pictures of the node numbers and contour plots of the solution you get at the moment. If your basis functions are correct, then test the rest of the code until the stiffness matrix symmetric matrix. Aug 18 '21 at 13:52
• No, the OP is correct in his assumption of the rigid body mode corresponding to a rotation about the axis between nodes 6 and 7. The correct stiffness matrix corresponding to this constraint set is singular; i.e. it has a single zero-eigenvalue with eigenvector representing the rigid body rotation. Aug 18 '21 at 15:02
• @BillGreene, you are correct! I tested thoroughly. The stiffness matrix has one zero eigenvalue even after fixing the nodes on of the edges. I learnt something new today. Aug 19 '21 at 22:10