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

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.

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

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.')
``````

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!

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.

You should collapse the subfunctions