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