Wrong results when using Hexahedron element

1

Hi
I am solving a nonlinear problem on a cube (Single element). Here is my minimal work:

from dolfin import *

n = 1
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), n,n,n)
#mesh = UnitCubeMesh.create(n,n,n,CellType.Type.hexahedron)

# Define boundary
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

#element for pressure field
Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)
#element for displacement field
Element2 = VectorElement("CG", mesh.ufl_cell(), 2)

# Defining the mixed function space
W_elem = MixedElement([Element1, Element2])
W = FunctionSpace(mesh, W_elem)
left = Left()
right = Right()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())

boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

dx = Measure('dx', domain=mesh, subdomain_data=subdomains, metadata={'quadrature_degree': 2})
ds = Measure('ds', domain=mesh, subdomain_data=boundaries, metadata={'quadrature_degree': 2})

n = FacetNormal(mesh)

# Fixed boundary condition at the left face
u_0 = Constant((0.0,0.0,0.0))
bc_left = DirichletBC(W.sub(1), u_0, boundaries, 1)

dw = TrialFunction(W)
v  = TestFunction(W)
w = Function(W)
p,u = split(w)

d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F

Ic = tr(C)
J  = det(F)

force  = Constant((1000,  0.0, 0.0))
# Total potential energy
psi = 0.4e6*(Ic - 3.)*dx - p*(J - 1)*dx - dot(force, u)*ds(2)
F1 = derivative(psi, w, TestFunction(W))

Jacob = derivative(F1, w, TrialFunction(W))
problem = NonlinearVariationalProblem(F1, w, bc_left, Jacob)

solver = NonlinearVariationalSolver(problem)
solver.solve()

(p, u) = w.split(True)

file = File("displacement.pvd")
file << u

The question is when I change the element type to hexahedron by using:

mesh = UnitCubeMesh.create(n,n,n,CellType.Type.hexahedron)

I get strange results that seems wrong. Here is a figure showing how the results are different:

11
11777×331 67.7 KB

Any idea what is wrong with the hexahedron element?
Thanks!

2

Removing your quadrature option from your integrals reduces the difference between the methods significantly.
Additionally, there is very limited support for visualizing hexahedrons (and especially higher order function spaces) in paraview.
You can compare you solutions (after removing the reduction of quadrature degree) by adding the following snippet:

for i,coord in enumerate(W.sub(0).collapse().tabulate_dof_coordinates()):
    print(coord, p.vector().get_local()[i])
print(p.vector().norm("l2"), u.vector().norm("l2"))

For the tetrahedon mesh:

[ 0.  0.  1.] 799361.270236
[ 0.  0.  0.] 799445.235254
[ 1.  1.  1.] 799879.027335
[ 1.  0.  1.] 799763.559526
[ 0.  1.  1.] 799024.881425
[ 1.  0.  0.] 799626.567357
[ 1.  1.  0.] 799763.559526
[ 0.  1.  0.] 799361.270236
2261407.290582405 0.0012195794459355432

For the hexahedral mesh:

[ 1.  0.  0.] 799720.082157
[ 0.  0.  0.] 799521.125508
[ 0.  1.  0.] 799521.125508
[ 1.  1.  0.] 799720.082157
[ 0.  0.  1.] 799521.125508
[ 1.  0.  1.] 799720.082157
[ 0.  1.  1.] 799521.125508
[ 1.  1.  1.] 799720.082157
2261668.6228875914 0.0012240027253240522

A similar comparison can be done for the displacements.

1 Like
3

Thanks. Is there any direct way to extract the nodal value at a specific coordinate when using Hexahedron element . For example I want to extract the nodal value by doing:

print(u(1,1,1))

But it returns:

*** Error:   Unable to intersect cell and point.
*** Reason:  Intersection is only implemented for simplex meshes.

It was an old issue and I thought it was fixed in 2019 version but it seems like it is still there.