2D result doesn't match from 3D code

General
1

Hi,
I am doing topology optimization.
I am trying to verify my 3D code (FEM) if it is working properly.
I am comparing strain energy from 2D code and 3D code (1 element along the z-axis).
The results seem very different.
I may be making a simple mistake.
Simplified code is as follows

##############3D################
from dolfin import *
E = Constant(1)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
Q = as_matrix([[lmbda+2*mu, lmbda, lmbda, 0, 0, 0],
                [lmbda, lmbda+2*mu, lmbda, 0, 0, 0],
                [lmbda, lmbda, lmbda+2*mu, 0, 0, 0],
                [0, 0, 0, mu, 0, 0],
                [0, 0, 0, 0, mu, 0],
                [0, 0, 0, 0, 0, mu],
              ])
length = 75 # [mm]
thickness =25 # [mm]
width = 1 # [mm]
resolution = 1 
f = Constant((0, -1, 0.0)) 
resX = int(resolution * length) # Num nodes in x axis
resY = int(resolution * thickness) # Num nodes in y axis
resZ = int(resolution * width) # Num nodes in z axis
Lx      = length
Ly      = thickness
Lz      = width
xmin    = 0.0
ymin    = 0.0
zmin    = 0.0
pmin    = Point(float(xmin), float(ymin), float(zmin))
xmax    = Lx
ymax    = Ly
zmax    = Lz
pmax    = Point(float(xmax), float(ymax), float(zmax))
Nx      = int(resX)
Ny      = int(resY)
Nz      = int(resZ)

mesh = BoxMesh.create([pmin, pmax], [resX, resY, resZ],CellType.Type.hexahedron)

V = VectorFunctionSpace(mesh, "CG", 1)

def strain_to_voigt(eps):
    return as_vector([eps[0, 0], eps[1, 1], eps[2, 2], 2*eps[0, 1], 2*eps[0, 2], 2*eps[1, 2] ])
def stress_from_voigt(S):
    ss = [[S[0], S[3], S[4]],
          [S[3], S[1], S[5]],
          [S[4], S[5], S[2]]]
    return as_matrix(ss)

def stress(u):
    return stress_from_voigt(Q*strain_to_voigt(sym(grad(u))))

def Stress(u):
    return stress_from_voigt(strain_to_voigt(sym(grad(u))))

def strain(u):
    return sym(grad(u))
class Symmetry(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < 1e-10

class BottomPoint(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] < 1e-10 and x[0] > 75-1e-10

class Loading(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[1] > 25-1e-10 and x[0] < 1.5

facets = MeshFunction('size_t', mesh, 1)
facets.set_all(0)
bottom = Loading()
bottom.mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)
u = TrialFunction(V)
v = TestFunction(V)
bc1 = DirichletBC(V.sub(0), Constant(0), Symmetry())
bc2 = DirichletBC(V.sub(1), Constant(0), BottomPoint(), method='pointwise')
bcs = [bc1, bc2]
u = TestFunction(V)
v = TrialFunction(V)
us = Function(V,name="Displacement")
dx = dx(metadata={'quadrature_degree': 3})
a = inner((stress(v)),strain(u))*dx
l =  dot(f, u)*ds(1)#+ dot(f, u)*dx 
solve(a == l,us,bcs=bcs , solver_parameters={"linear_solver" : "lu", "preconditioner" : "gamg"})
cost = assemble(inner((stress(us)), strain(us))*dx) #strain energy

#####################2D##################################
E = Constant(1)
nu = Constant(0.3)

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
Q = as_matrix([[lmbda+2*mu, lmbda, 0],
                [lmbda, lmbda+2*mu, 0],
                [0, 0, mu]])
def strain_to_voigt(e):
    return as_vector((e[0, 0], e[1, 1], 2*e[0, 1]))

def stress_from_voigt(sigma_voigt):
    return as_matrix(((sigma_voigt[0], sigma_voigt[2]), (sigma_voigt[2], sigma_voigt[1])))

def stress(u):
    return stress_from_voigt(Q*strain_to_voigt(sym(grad(u))))

def strain(u):
    return sym(grad(u))
    
class Symmetry(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < 1e-10

class BottomPoint(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] < 1e-10 and x[0] > 75-1e-10

class Loading(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[1] > 25-1e-10 and x[0] < 1.5
    
length = 75
width = 25
res = 1

Lx      = length
Ly      = width
xmin    = 0.0
ymin    = 0.0
pmin    = Point(float(xmin), float(ymin))
xmax    = Lx
ymax    = Ly
pmax    = Point(float(xmax), float(ymax))
resX = int(res * length) # Num nodes in x axis
resY = int(res * width) # Num nodes in y axis
Nx      = int(resX)
Ny      = int(resY)

mesh = RectangleMesh.create([pmin, pmax], [resX, resY],CellType.Type.quadrilateral)





X = FunctionSpace(mesh, "CG", 1)
V = VectorFunctionSpace(mesh, "CG", 1)
N = Function(X).vector().size()
problemSize = N*3

facets = MeshFunction('size_t', mesh, 1)
facets.set_all(0)
bottom = Loading()
bottom.mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)
bc1 = DirichletBC(V.sub(0), Constant(0), Symmetry())
bc2 = DirichletBC(V.sub(1), Constant(0), BottomPoint(), method='pointwise')

bcs = [bc1, bc2]
  
#     return bcs, ds
u = TrialFunction(V)
v = TestFunction(V)
us = Function(V,name="Displacement")
a = inner((stress(u)),strain(v))*dx
L = inner(f, v)*ds(1)
solve(a == L,us,bcs=bcs , solver_parameters={"linear_solver" : "lu", "preconditioner" : "gamg"})
cost = assemble(inner((stress(us)), strain(us))*dx)

Sorry for the long code. For 3D it’s 771 vs 2D 120