Convergence with CG, 1 elements for incompressible elasticity

dolfinx
1

Hi all,

I adapted a code from a previous forum post to the latest version of fenicsx. the simulation runs smoothly for 2nd order continuous galerkin elements, but fails to converge for first order elements, how is that so? this was done using
dolfinx version 0.8

import numpy as np
import ufl

from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot, io, nls, log, cpp
import basix
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

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

P2 = basix.ufl.element(
    "CG", 
    domain.topology.cell_name(), 
    1, 
    shape=(domain.geometry.dim, )
)

P1 = basix.ufl.element(
    "CG", 
    domain.topology.cell_name(), 
    1,
)

# combine and define our mixed element
mixed = basix.ufl.mixed_element([P2, P1])
state_space = fem.functionspace(domain, mixed)

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 = NonlinearProblem(R, state, bcs)
solver = 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()

in particular, changing the displacement element to first order:

P2 = basix.ufl.element(
    "CG", 
    domain.topology.cell_name(), 
    1, 
    shape=(domain.geometry.dim, )
)

then fails the simulation, as such:

RuntimeError: Newton solver did not converge because maximum number of iterations reached

thanks in advance.

2

With incompressibility constraints you are not free to arbitrarily choose the degrees of the spaces in your mixed space. For incompressible fluid dynamics, a keyword to search for is inf-sup stability; incompressible elasticity has a similar concept.

Knowledge of these concepts should come from any FEM graduate course: in this forum we can’t replace a course, just point you in the right direction.