Fully incompressible formulation for a non-linear hyperelasticity problem in FenicsX

dolfinx
1

Hi everyone,

I am trying to convert a code written in Fenics to FenicsX. Here we use a fully incompressible formulation, solving for the displacement and the hydrostatic pressure simultaneously. However, it doesn’t work (aka, doesn’t converge), and I’m not able to see why.

I’ve tried to adapt this code, which runs fine as it is, to my case: Hyperelasticity — FEniCSx tutorial

The code reads as follows:

import numpy as np
import ufl

from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot, io, nls, log, cpp

domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [1, 1, 1]], [2, 2, 2], mesh.CellType.hexahedron)

P2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)

state_space = fem.FunctionSpace(domain, P2 * P1)
V, _ = state_space.sub(0).collapse()

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], 1)

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

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
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=PETSc.ScalarType)

left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

T = fem.Constant(domain, PETSc.ScalarType((0, 0, -1.5)))

state = fem.Function(state_space, name="state")
test_state = ufl.TestFunctions(state_space)

u, p = state.split()
v, q = test_state

# Kinematics
d = len(u)
I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))

# Elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.5)
mu = fem.Constant(domain, E/(2*(1 + nu)))

# Stored strain energy density (fully incompressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J)

# Hyper-elasticity
P = ufl.diff(psi, F) + p * J * ufl.inv(F.T)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)

R = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, T)*ds(2) + q * (J - 1) * dx

problem = fem.petsc.NonlinearProblem(R, state, bcs)
solver = nls.petsc.NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

solver.solve(state)

If I implement the same in Fenics it does run nicely, so I think my error is in the FenicsX syntax. Is anyone able to see what could cause the non-convergent behavior?

2

I think the error is in your Dirichlet bc. You should not use the collapsed subspace, i.e.

is wrong. Instead, you should use the non-collapsed subspace for Dirichlet bc creation, as described in my previous answer, e.g.:

left_dofs = fem.locate_dofs_topological((state_space.sub(0), V), facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, state_space.sub(0))]

In your current code, because you use the collapsed subspace, you are effectively trying to solve the problem without any Dirichlet bcs, which would understandable give rise to the convergence issue.

1 Like
3

Hi @conpierce8, and thanks a lot for answering!

That makes a lot of sense, thanks for the clarification. I completely get why it won’t work. However, if I try replacing the two lines with yours, I get an error message I don’t understand …

$ python mwe_pressure.py 
Traceback (most recent call last):
  File "/home/aashild/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/bcs.py", line 125, in __init__
    super().__init__(_value, dofs, V)  # type: ignore
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: numpy.ndarray[numpy.int32])
    4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: List[numpy.ndarray[numpy.int32][2]], V: dolfinx::fem::FunctionSpace)

Invoked with: array([0., 0., 0.]), [array([  0,   1,   2,   6,   7,   8,  12,  13,  14,  18,  19,  20,  27,
        28,  29,  30,  31,  32,  42,  43,  44,  51,  52,  53,  66,  67,
        68, 147, 148, 149, 153, 154, 155, 159, 160, 161, 168, 169, 170,
       174, 175, 176, 186, 187, 188, 205, 206, 207, 211, 212, 213, 217,
       218, 219, 223, 224, 225, 232, 233, 234, 244, 245, 246, 339, 340,
       341, 345, 346, 347, 351, 352, 353, 360, 361, 362], dtype=int32), array([  0,   1,   2,  36,  37,  38,  45,  46,  47,  51,  52,  53,   9,
        10,  11,  12,  13,  14,  39,  40,  41,  33,  34,  35,  15,  16,
        17,  90,  91,  92, 114, 115, 116,  93,  94,  95,  99, 100, 101,
       117, 118, 119, 105, 106, 107,  54,  55,  56,  78,  79,  80,  57,
        58,  59,  81,  82,  83,  63,  64,  65,  69,  70,  71, 126, 127,
       128, 129, 130, 131, 132, 133, 134, 138, 139, 140], dtype=int32)], FunctionSpace(Mesh(VectorElement(FiniteElement('Q', hexahedron, 1), dim=3), 0), VectorElement(FiniteElement('Q', hexahedron, 2), dim=3))

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/aashild/Documents/Projects/EMI/mwe_pressure.py", line 35, in <module>
    bcs = [fem.dirichletbc(u_bc, left_dofs, state_space.sub(0))]
  File "/home/aashild/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/bcs.py", line 182, in dirichletbc
    return formcls(value, dofs, V)
  File "/home/aashild/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/bcs.py", line 127, in __init__
    super().__init__(_value, dofs, V._cpp_object)  # type: ignore
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.DirichletBC_float64(g: numpy.ndarray[numpy.float64], dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    2. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Constant<double>, dofs: numpy.ndarray[numpy.int32], V: dolfinx::fem::FunctionSpace)
    3. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: numpy.ndarray[numpy.int32])
    4. dolfinx.cpp.fem.DirichletBC_float64(g: dolfinx::fem::Function<double>, dofs: List[numpy.ndarray[numpy.int32][2]], V: dolfinx::fem::FunctionSpace)

Invoked with: array([0., 0., 0.]), [array([  0,   1,   2,   6,   7,   8,  12,  13,  14,  18,  19,  20,  27,
        28,  29,  30,  31,  32,  42,  43,  44,  51,  52,  53,  66,  67,
        68, 147, 148, 149, 153, 154, 155, 159, 160, 161, 168, 169, 170,
       174, 175, 176, 186, 187, 188, 205, 206, 207, 211, 212, 213, 217,
       218, 219, 223, 224, 225, 232, 233, 234, 244, 245, 246, 339, 340,
       341, 345, 346, 347, 351, 352, 353, 360, 361, 362], dtype=int32), array([  0,   1,   2,  36,  37,  38,  45,  46,  47,  51,  52,  53,   9,
        10,  11,  12,  13,  14,  39,  40,  41,  33,  34,  35,  15,  16,
        17,  90,  91,  92, 114, 115, 116,  93,  94,  95,  99, 100, 101,
       117, 118, 119, 105, 106, 107,  54,  55,  56,  78,  79,  80,  57,
        58,  59,  81,  82,  83,  63,  64,  65,  69,  70,  71, 126, 127,
       128, 129, 130, 131, 132, 133, 134, 138, 139, 140], dtype=int32)], <dolfinx.cpp.fem.FunctionSpace object at 0x7f26038043f0>

I’ve tried variants of this (e.g., passing in left_dofs[0] instead of left_dofs) but so far without any success. Does the above error message by any chance make more sense to you?

4

My mistake. I should have tried running my proposed solution before recommending it. Yes, the error message makes sense. When you call locate_dofs_topological with a list of two function spaces, it returns a list of two ndarrays holding the dof indices for each function space. fem.dirichletbc is therefore called with a list of two ndarrays for its second argument, which is usage 4 in the error message. For that usage, the only allowable type for the first argument is a fem.Function (Constants and ndarrays are not allowed). Therefore, what you need to do is:

  1. Create an auxiliary Function on the function space V. This function will hold the value you want to prescribe in your bc.
  2. Interpolate the constant value (in your case, (0.0, 0.0, 0.0)) into the function you just created.
  3. Pass the Function to dirichletbc as the first argument.

When you make these changes, you’ll probably find that you still have nonconvergence. I know the solution to this as well, but I don’t know why it’s the solution. :confused: In general, you need to use ufl.split(state) when you create functions/testfunctions for use in forms, and state.split() when you split your solution into components for post-processing. Here’s a code that converges and appears to produce the correct result:

import numpy as np
import ufl

from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot, io, nls, log, cpp

domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [1, 1, 1]], [2, 2, 2], mesh.CellType.hexahedron)

P2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)

state_space = fem.FunctionSpace(domain, P2 * P1)
V, _ = state_space.sub(0).collapse()

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], 1)

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

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

u0_bc = fem.Function(V)
u0 = lambda x: np.zeros_like(x, dtype=PETSc.ScalarType)
u0_bc.interpolate(u0)

left_dofs = fem.locate_dofs_topological((state_space.sub(0), V), facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u0_bc, left_dofs, state_space.sub(0))]

T = fem.Constant(domain, PETSc.ScalarType((0, 0, -1.5)))

state = fem.Function(state_space, name="state")
test_state = ufl.TestFunctions(state_space)

u, p = ufl.split(state)
v, q = test_state

# Kinematics
d = len(u)
I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))

# Elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.5)
mu = fem.Constant(domain, E/(2*(1 + nu)))

# Stored strain energy density (fully incompressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J)

# Hyper-elasticity
P = ufl.diff(psi, F) + p * J * ufl.inv(F.T)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)

R = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, T)*ds(2) + q * (J - 1) * dx

problem = fem.petsc.NonlinearProblem(R, state, bcs)
solver = nls.petsc.NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

solver.solve(state)
u_sol, p_sol = state.split()
2 Likes
5

Thank you so much @conpierce8 - worked like a charm! I’m totally ok with using ufl.split without understanding why :slight_smile:

Also, the link to your other post will be very useful – it is actually my next step to do this component-wise!

1 Like
6

Hello,

I tried add one more displacement load on the right face in the follow, and make the traction and body force as 0.

### Add displacement on u_x right 
U_right_x = 1
facet_right = mesh.locate_entities_boundary(domain, fdim, right)
dofs_right_x = fem.locate_dofs_topological(V.sub(1), domain.topology.dim-1, facet_right)
bc_right = fem.dirichletbc(ScalarType(U_right_x), dofs_right_x, V.sub(1))

bcs = [bc_left, bc_right]

However, it doesn’t converge with the error “Newton solver did not converge because maximum number of iterations reached”. Do you know why?

7

Could you post your complete code? Without it, I would be completely unable to say why it doesn’t converge.

8

Thanks for your help in advance, here it is:

import numpy as np
import ufl

from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot, io, nls, log, cpp
from petsc4py.PETSc import ScalarType

domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [1, 1, 1]], [1, 1, 1], mesh.CellType.hexahedron)

P2 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 1)
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)

state_space = fem.FunctionSpace(domain, P2 * P1)
V, _ = state_space.sub(0).collapse()

def left(x):
    return np.isclose(x[0], 0)

def right(x):
    return np.isclose(x[0], 1)

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

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

### original boundary condition
u1_bc = fem.Function(V)
u1 = lambda x: np.zeros_like(x, dtype=PETSc.ScalarType)
u1_bc.interpolate(u1)
left_dofs = fem.locate_dofs_topological((state_space.sub(0), V), facet_tag.dim, facet_tag.find(1))
bc_left = fem.dirichletbc(u1_bc, left_dofs, state_space.sub(0))
bcs = [bc_left]

The following is the displacement u_x applied on the right face:

### Add displacement u_x on right face

u_x = ScalarType(1)
facet_right = mesh.locate_entities_boundary(domain, fdim, right)
dofs_right_x = fem.locate_dofs_topological((V.sub(0),V), domain.topology.dim-1, facet_right)
bc_right = fem.dirichletbc(ux, dofs_right_x, V.sub(0))

bcs = [bc_left, bc_right]

The last part is the same with the original code, I just set the traction on the boundary as 0 here.

T = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))

state = fem.Function(state_space, name="state")
test_state = ufl.TestFunctions(state_space)

u, p = ufl.split(state)
v, q = test_state

# Kinematics
d = len(u)
I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J  = ufl.variable(ufl.det(F))

# Elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.5)
mu = fem.Constant(domain, E/(2*(1 + nu)))

# Stored strain energy density (fully incompressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J)

# Hyper-elasticity
P = ufl.diff(psi, F) + p * J * ufl.inv(F.T)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)

R = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, T)*ds(2) + q * (J - 1) * dx

problem = fem.petsc.NonlinearProblem(R, state, bcs)
solver = nls.petsc.NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

solver.solve(state)
u_sol, p_sol = state.split()
9

You should not use V.sub(0), but rather state_space.sub(0).sub(0) as I explained above.

Since you are applying a displacement on the right boundary, you should not apply a traction there, i.e. you need to remove - ufl.inner(v, T)*ds(2) from the weak form R.

2 Likes
10

Thanks a lot for your explanation.

Could I ask one more basic question here? I am confused with the way of setting dofs, is there any difference between

fem.locate_dofs_topological()

and

fem.locate_dofs_geometrical()

as I thought they are the same to locate dofs.

11

See Difference between locate_dofs_topological and locate_dofs_geometrically - #2 by dokken

1 Like