Hyper elasticity tutorial in 2D and 3D comparison

General
1

Hello everyone, I am trying to run the hyper-tutorial example Hyperelasticity — FEniCSx tutorial in cases of 2D and 3D with modification only in the neo-Hookean energy model. I tried to run for case of no external load condition and tested. In 3D it works perfectly fine and produces no displacement but in 2D with the same setup it produces strange displacement values.
The code with both cases of 2D and commented 3D is shown below,

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl

from mpi4py import MPI
from dolfinx import fem, mesh, plot
L = 20
                                #############For 3D case#######################
# domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)
# B = fem.Constant(domain, default_scalar_type((0, 0,0)))
# T = fem.Constant(domain, default_scalar_type((0, 0,0)))

                                #############For 2D case#######################
domain = mesh.create_rectangle (comm=MPI.COMM_WORLD, points=((0.0, 0.0), (20.0,1.0)), n=(1,1), cell_type=mesh.CellType.quadrilateral)
B = fem.Constant(domain, default_scalar_type((0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0)))

element = ufl.VectorElement("Lagrange", domain.ufl_cell(),1)
V=fem.FunctionSpace(domain, element)
def left(x):
    return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets])
marked_values = np.hstack([np.full_like(left_facets, 1)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))

bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

v = ufl.TestFunction(V)
u = fem.Function(V)

# Spatial dimension
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = ((mu/2) * (((J**(-2./3.)) * Ic) - 3)) + ((lmbda/2) * (((J**2.) /2) - (1/2) - ufl.ln(J)) )
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds

problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
solver.solve(u)
print (u.x.array)

The boundary condition for both the case is almost the same, major change is in the dimensions in forces. Could anyone please confirm me where it is going wrong, as I am unable to figure where it is going wrong for past one week. Thanks everyone for your support and time.

Any tips on the above topic is greatly appreciated. Thanks again for your considerations.

2

If you want anyone to help you, you have to be more precise in what behavior you see and what you would have expected instead.

domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)
domain = mesh.create_rectangle (comm=MPI.COMM_WORLD, points=((0.0, 0.0), (20.0,1.0)), n=(1,1), cell_type=mesh.CellType.quadrilateral)

To double check: are you aware that you are using a mesh with a single cell? If so, can you explain why is it that you expect to have an accurate solution with such a coarse mesh?

3

Thanks for your response Mr.Ballarin. Primarily I am checking here just the behavior of the problem in no load condition.
I am expecting that the Displacement have to be zero, but unfortunately in 2D case, without any changes the solver shows values for u. Ideally it should be zero in no load condition which is perfectly working in 3D case.

So I am confused which part of the code is actually going wrong.

Thanks for your support and consideration.

4

IMHO, the error must be in

psi = ((mu/2) * (((J**(-2./3.)) * Ic) - 3)) + ((lmbda/2) * (((J**2.) /2) - (1/2) - ufl.ln(J)) )

in a way that when you do

P = ufl.diff(psi, F)

the result is the sum between a term that depends on u and a term which does not depend on u, and the latter ends up as a forcing term to your problem.

Unfortunately hyperelasticity is not my field of research, so I can’t be of much more help. Are you sure that the expression for P is valid both in 2D and 3D? Is it possible that maybe some of the coefficients (like the 2/3 exponent) are slightly different in 2D?

5

No Mr.Ballarin, the psi expression for the Neo-Hookean energy model is not different for both the cases.
Both the term having the expression J that depends on u. Is that the psi expression you were mentioning in the question?
Thanks for your time and considerations.

6

OK, then you’ll have to wait for someone who is more expert in the field. My impression is that there is a forcing term coming in through P (and it shouldn’t), but I can’t be sure if that is true, or pinpoint why it happens.

7

The problem comes from the fact that the free energy expression is valid only for the 3D case. As suggested by @francesco-ballarin, with your expression and a 2D gradient you end up with forcing terms. To avoid this, you can just change your script by using 3D terms in a 2D formulation such as:

# Identity tensor
I = ufl.variable(ufl.Identity(3))
# Deformation gradient
def grad_3D(u):
    return ufl.as_matrix([[u[0].dx(0), u[0].dx(1), 0], 
                          [u[1].dx(0), u[1].dx(1), 0], 
                          [0, 0, 0]])
F = ufl.variable(I + grad_3D(u))
8

Thanks a lot for your response Mr. Bleyerj. But unfortunately this resulting in the shape mismatch error in the weak form equation. Changing here as 3D for the deformation gradient makes the mismatch in the weak form when solving in 2D case(shape mismatch between P and grad(v) for the inner operation). Additionally any terms required these for these adaptation. Could you please confirm me here?. Thanks again for your time and considerations.
The error message after the above adaptation is

Traceback (most recent call last):
  File "/mnt/e/PPP/Project/Chapter 1-29-11/3 point bending setup/hyper.py", line 81, in <module>
    F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds
  File "/usr/lib/python3/dist-packages/ufl/operators.py", line 167, in inner
    return Inner(a, b)
  File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 162, in __new__
    raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <Grad id=140324846198336> and <VariableDerivative id=140325375135232>
9

Obviously you have to use grad_3D instead of ufl.grad in the weak form as well

1 Like
10

Hello. I have a small question in visualising the value of deformation gradient, I use a method like creating again a function space for visualising it as shown below in code. I am not sure whether any other approaches available. But here since the F given as 3X3 for 2D case, while visualising with this method yields an error for shape mismatch.

#Printing the deformation gradient
element2 = ufl.TensorElement("DG",domain.ufl_cell(),0)
V2=fem.FunctionSpace(domain,element2)
s = ufl.TrialFunction(V2)
t = ufl.TestFunction(V2)

a = ufl.inner(s, t)*dx
L = ufl.inner(F, t)*dx #Here F is the deformation gradient 

problem2 = LinearProblem(a,L)
deformation = problem2.solve()

val=deformation.x.array
print (val)

So am I making anything wrong or any other better approach for visualising the deformation gradient. Could anyone please tell me about the above regard. Thanks for your responses.

11

Set the dim to (3,3) as a keyword argument.