Mixed (u/p) hyperelasticity: Problem with component-wise BCs

dolfinx
1

Dear all,

I am trying to solve the equations of incompressible hyperelasticity in dolfinx 0.7.3 by means of a mixed (hybrid u/p) element formulation. In doing so, I am facing problems with the prescription of the BCs. More precisely, I try to prescribe the x-component of the displacement vector on the right edge of the domain. However, when visualising the results in paraview, it turns out that not only u_x has been set to this value, yet all three components.In other words, the following code snippet doesn’t seem to work appropriately; please find an MWE below.

state_space = fem.FunctionSpace(domain, P2 * P1)
 # ...
right_u_dofs = fem.locate_dofs_geometrical((state_space.sub(0).sub(0),u_x_space), right)
bc_right = fem.dirichletbc(u_x_traction, right_u_dofs, state_space.sub(0).sub(0))

The problem is very similar to recent posts (Incompressible Mooney-Rivlin with Component Wise BCs and Fully incompressible formulation for a non-linear hyperelasticity problem in FenicsX - #5 by aashildte). However, unfortunately, the solutions that have been posted there do not seem to work with my version of fenicsx.

Thanks in advance for your help!

import numpy as np
import ufl

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

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 = ufl.VectorElement("CG", domain.ufl_cell(), 2)
P1 = ufl.FiniteElement("CG", domain.ufl_cell(), 1)

state_space = fem.FunctionSpace(domain, P2 * P1)
u_space, _ = 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

# Fix displacement DOF on left face
u0_bc = fem.Function(u_space)
u0 = lambda x: np.stack((np.zeros(x.shape[1]), np.zeros(x.shape[1]), np.zeros(x.shape[1])))
u0_bc.interpolate(u0)

left_dofs = fem.locate_dofs_geometrical((state_space.sub(0),u_space), left)
leftBc = fem.dirichletbc(u0_bc, left_dofs, state_space.sub(0))

# non-zero x displacement on right face -- there seems to be a problem, here
u_x_space, subsubmap = u_space.sub(0).collapse()
u_x_traction = fem.Function(u_x_space)

def incremented_displacement_expression(x):
    return np.full(x.shape[1], 1.0e-1)
u_x_traction.interpolate(incremented_displacement_expression)

right_u_dofs = fem.locate_dofs_geometrical((state_space.sub(0).sub(0),u_x_space), right)
bc_right = fem.dirichletbc(u_x_traction, right_u_dofs, state_space.sub(0).sub(0))

bcs = [leftBc,bc_right]

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}
dx = ufl.Measure("dx", domain=domain, metadata=metadata)

Rp =  q * (J - 1) * dx 

R = ufl.inner(ufl.grad(v), P)*dx  + Rp


log.set_log_level(log.LogLevel.INFO)

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()

t = 1.
### Incremental post-processing ###
print('Performing paraView data export...')
u_sol, p_sol = state.split() #u,JStar,pStar = ufl.split(up)
# paraview export
u_sol.name = "Displacement"
p_sol.name = "Pressure"
# paraview data
with io.VTXWriter(domain.comm, 'resultsMWE_discGr/disp_inc.bp', [u_sol.sub(0),u_sol.sub(1),u_sol.sub(2)], engine="BP4") as vtx: # write displacement solution
	vtx.write(t)
with io.VTXWriter(domain.comm, 'resultsMWE_discGr/pStar_inc.bp', [p_sol], engine="BP4") as vtx: # pStar
	vtx.write(t)
print('paraView export has been done.')
2

Update: I also tried this in Fenicsx 0.8.0, with the definition of the “mixed” function space changed.

P2 = uflx.element("Lagrange", domain.basix_cell(), 2, shape=(domain.geometry.dim,))
P1 = uflx.element("Lagrange", domain.basix_cell(), 1)

TH = uflx.mixed_element([P2, P1])
state_space = fem.functionspace(domain, TH)

However, the problem remains the same. I would highly appreciate any help for correctly defining Dirichlet BCs for components of the vector-valued function! :confused:

3

When trying to fix the issues, I found out that this is not a problem of the definition of the BCs. In this regards, everything is working fine. However, there actually is a problem with the Paraview VTX output here, into which the values for the different displacement components do not have been written correctly.

4

You should collapse the subfunctions