# Displacement boundary condition issue for compression test in mixed formulation

Hello everyone,

To explain the scenario, first I used a hyperelastic model like the one from the tutorials to perform a compression test simulation. It is being performed on a 3D cylinder geometry. The displacement of 0.001m/s is applied for 3 seconds in z direction and then it is kept constant for a total time of 1000 seconds.

However, adapting the same boundary conditions is giving me errors and I am stuck with it for a few days now. Could anyone help me understand how to adapt it for mixed formulation or what am I doing wrong?

So using the hyperelastic model from tutorials which works fine. I used the following boundary conditions:

def bc_vec_zero(x):
return np.zeros((dim, x.shape[1]), dtype=PETSc.ScalarType)

list_bc = []

u_bc = np.array((0,) * dim, dtype=default_scalar_type)

# Fix bottom surface with facet_tag = 3
bottom_facets = facet_tag.indices[facet_tag.values == 3]
bottom_dofs = fem.locate_dofs_topological(V, fdim, bottom_facets)
list_bc.append(fem.dirichletbc(u_bc, bottom_dofs, V))

#Strain function for top surface
axial_displacement = 0.001

def axial_strain(t):
if t <= 3:
return -t * axial_strain_rate
else:
return -0.003
# Assign boundary condition to top surface with facet tag = 2
u_bc_2 = fem. Constant(domain, ScalarType(axial_strain(t)))
top_facets = facet_tag.indices[facet_tag.values == 2]
top_dofs = fem.locate_dofs_topological(V.sub(2), fdim, top_facets)
list_bc.append(fem.dirichletbc(u_bc_2, top_dofs, V.sub(2)))


Now I want to use the same test with mixed formulation. Although I am not able to solve this issue for a few days now. This is how I am trying to adapt to the displacement boundary conditions:

#-----------------------Function Spaces and Mixed Space-------------------#

u_el = ufl.VectorElement("CG", domain.ufl_cell(), 2) # u quad element
p_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # p linear element
element = ufl.MixedElement([u_el, p_el]) # Mixed Element

V = fem.FunctionSpace(domain, element) # Mixed Function space for all dofs
up = fem.Function(V) # total dofs in the problem/ func for all dofs
print("up dofs:", len(up.x.array))

u_sp = fem.FunctionSpace(domain, u_el) # u Function space
p_sp = fem.FunctionSpace(domain, p_el) # p Function space

Vu, up_to_u = V.sub(0).collapse() # dofs in u func subspace
u_n = fem.Function(Vu) # dofs in u func
#print("u_n dofs:", len(u_n.x.array))

Vp, up_to_p = V.sub(1).collapse() # dofs in u func subspace
p_n = fem.Function(Vp) # dofs in u func
#print("p_n dofs:", len(p_n.x.array))

u, p = ufl.split(up)
print("all dofs:", len(u))
(del_u, del_p) = ufl.TestFunctions(V)
#print("all dofs:", len(del_u))

def bc_vec_zero(x):
return np.zeros((dim, x.shape[1]), dtype=PETSc.ScalarType)
list_bc = []

# Displacement Boundary
u_bc = fem.Function(u_sp)
u_bc.interpolate(bc_vec_zero)

#-----------------Fixed Bottom (3)---------------------

facets = facet_tag.indices[facet_tag.values == 3] # fixed in x
left_dofs = fem.locate_dofs_topological((V.sub(0).sub(0), u_sp.sub(0)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, left_dofs, V.sub(0)))

facets = facet_tag.indices[facet_tag.values == 3] # fixed in y
bottom_dofs = fem.locate_dofs_topological((V.sub(0).sub(1), u_sp.sub(1)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, bottom_dofs, V.sub(0)))

facets = facet_tag.indices[facet_tag.values == 3] # fixed in z
bottom_dofs = fem.locate_dofs_topological((V.sub(0).sub(2), u_sp.sub(2)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, bottom_dofs, V.sub(0)))

#--------- Compressive Strain on Top Surface (2)------------#

displacement_rate = 0.001 # 1mm/s per second

def axial_strain(t):
if t <= 3.0:
return -t * displacement_rate
else:
return -0.003

u_bc_2 = fem. Constant(domain, ScalarType(axial_strain(t)))
top_facets = facet_tag.indices[facet_tag.values == 2]
top_dofs = fem.locate_dofs_topological((V.sub(0).sub(2), u_sp.sub(2)), fdim, top_facets)
list_bc.append(fem.dirichletbc(u_bc_2, top_dofs, V.sub(0)))



The screenshot of the error is:

Without a reproducible example it is hard to give you guidance. The error message stems from:

Which is due to the fact that you are sending in the wrong dofs.
It should simply be:

u_bc_2 = fem. Constant(domain, ScalarType(axial_strain(t)))
top_facets = facet_tag.indices[facet_tag.values == 2]
top_dofs = fem.locate_dofs_topological(V.sub(0).sub(2), fdim, top_facets)
list_bc.append(fem.dirichletbc(u_bc_2, top_dofs, V.sub(0).sub(2)))


I am not 100 % sure that this will run, so you would have to get back to me with a minimal reproducible example if it doesn’t

Thank you Dokken for the solution. It works perfectly fine now. I have attached the code at the end for others. Although I have 2 questions.

1. I see that I did not understand this correctly. Could you shortly explain (or point me to the right source/documentation) why it has to be like you mentioned? I usually fix, eg: displacement in x direction in the following way. Therefore, I tried it in a similar way. Why does the following way did not work.
# bottom surface (facet tag 3) of cylinder in x direction
facets = facet_tag.indices[facet_tag.values == 3] # fixed in x
left_dofs = fem.locate_dofs_topological((V.sub(0).sub(0), u_sp.sub(0)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, left_dofs, V.sub(0)))

1. I want to calculate stress/force on the top surface of the cylinder so that I can analyse the stress strain relations in the end. The idea was to calculate the traction on the top surface and multiply it by the area of the circular surface. For this, I added the code mentioned below in the time loop and it calculates it. However I am facing issues to save this in a file. I usually use fem. Expression() to evaluate and save data using VTXWriter. But I see that Expressions won’t work for integral measures. What would be the right way to save this data in an external file?
    #---------stress on top surface-------#
traction = ufl.inner(df, ufl.dot(P,normal_vec))
Top_force = fem.assemble_scalar(fem.form(traction * ds(2)))
Top_stress = Top_force / (3.14 * pow(0.105,2))
print('traction =',Top_stress)


Here is the code that works fine now. Instead of NeoHookean, I have used a model called Yeoh to represent elasticity as I have the material parameters available for the compression test I am trying to simulate. It is a poroelastic problem.

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

import numpy as np
import ufl
import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

from mpi4py import MPI
from dolfinx import fem
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc #solvers
from pathlib import Path

from dolfinx.io import gmshio
import meshio
#------------------------------meshio-------------------------------#
folder = "yeoh_bi"
#---------------------- 2nd order tri elements (set order 1 in gmsh)-------------------#
proc = MPI.COMM_WORLD.rank

# Extract data and returns meshio mesh including physical markers
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]})
return out_mesh

if proc == 0:
print(msh)
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "tetra")
line_mesh = create_mesh(msh, "triangle")
meshio.write(Path(folder)/"mesh.xdmf", triangle_mesh)
meshio.write(Path(folder)/"mt.xdmf", line_mesh)
MPI.COMM_WORLD.barrier()

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, Path(folder)/"mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid") # cell tags
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, Path(folder)/"mt.xdmf", "r") as xdmf:
facet_tag = xdmf.read_meshtags(domain, name="Grid") # facet tags

fdim = domain.topology.dim - 1
dim = domain.topology.dim

#-----------------------Function Spaces and Mixed Space for displacement and pressure-------------------#

u_el = ufl.VectorElement("CG", domain.ufl_cell(), 2) # u quad element
p_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # p linear element
element = ufl.MixedElement([u_el, p_el]) # Mixed Element

V = fem.FunctionSpace(domain, element) # Mixed Function space for all dofs
up = fem.Function(V) # total dofs in the problem/ func for all dofs

u_sp = fem.FunctionSpace(domain, u_el) # u Function space
p_sp = fem.FunctionSpace(domain, p_el) # p Function space

Vu, up_to_u = V.sub(0).collapse() # dofs in u func subspace
u_n = fem.Function(Vu) # dofs in u func

Vp, up_to_p = V.sub(1).collapse() # dofs in u func subspace
p_n = fem.Function(Vp) # dofs in u func

u, p = ufl.split(up)
print("all dofs:", len(u))
(del_u, del_p) = ufl.TestFunctions(V)

#---------------------------temporal parameters---------------------------#
# Define temporal parameters
t = 0  # Start time
T = 100 #80000
num_steps = 100
dt = T / num_steps  # time step size

#---------------------------Initial Conditions---------------------------#

# Initialize history values
up.x.array[:] = 0.0
u_n.x.array[:] = 0.0

#---------------------------Boundary Conditions---------------------------#

def bc_vec_zero(x):
return np.zeros((dim, x.shape[1]), dtype=PETSc.ScalarType)

def bc_scalar_zero(x):
return np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)

list_bc = []

# Displacement Boundary
u_bc = fem.Function(u_sp)
u_bc.interpolate(bc_vec_zero)

#-----------------Fixed Bottom (3)---------------------

facets = facet_tag.indices[facet_tag.values == 3] # fixed in x
left_dofs = fem.locate_dofs_topological((V.sub(0).sub(0), u_sp.sub(0)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, left_dofs, V.sub(0)))

facets = facet_tag.indices[facet_tag.values == 3] # fixed in y
bottom_dofs = fem.locate_dofs_topological((V.sub(0).sub(1), u_sp.sub(1)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, bottom_dofs, V.sub(0)))

facets = facet_tag.indices[facet_tag.values == 3] # fixed in z
bottom_dofs = fem.locate_dofs_topological((V.sub(0).sub(2), u_sp.sub(2)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, bottom_dofs, V.sub(0)))

#--------- Compressive Strain on Top Surface (2)------------#

displacement_rate = 0.001 # 1mm/s per second

def axial_strain(t):
if t <= 3.0:
return -t * displacement_rate
else:
return -0.003

def u_exact(x):
return np.arry([0, 0, -0.001], (dim, x.shape[1]), dtype=PETSc.ScalarType)

u_bc_2 = fem. Constant(domain, ScalarType(axial_strain(t)))
top_facets = facet_tag.indices[facet_tag.values == 2]
top_dofs = fem.locate_dofs_topological(V.sub(0).sub(2), fdim, top_facets)
list_bc.append(fem.dirichletbc(u_bc_2, top_dofs, V.sub(0).sub(2)))

#----Normal Vector--------#

normal_vec = ufl.FacetNormal(domain)

#------------------------------Parameters-----------------------------------#

# Parameters for hyperelastic Yeoh model (could be neohookean as well)
c1 = 0.29 #material parameters
c2 = 1.01 #material parameter
d = 1.60 #material parameter
e= 0.90 #void ratio
k = 1.1e-10 #permeability

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
b = ufl.variable(F * F.T)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Time Derivatives
vel_s = (u - u_n) / dt # Solid Velocity

#Void ratio
nF = e/(1+e)

#--------------- Darcy Law --------------#
w_FS = dar/nF

#-----Yeoh stress (could be Neo Hookean as well)--------#

P_yeoh_b = ((c1) + (c2 * 2 * (Ic - 3))) * (1/J)
P_yeoh_i = ((1/d) * (J - 1))
P_yeoh = (P_yeoh_b * b) + (P_yeoh_i * I)

#------Total Stress----#
P = P_yeoh - p*I

# Integration measures

#-----Weak Forms--------#

#Momentum Balance
Mom_B = ufl.inner(P, ufl.grad(del_u)) * dx # Solid Stress part

#Volume Balance
Vol_M = (ufl.div(vel_s) * del_p) * dx # div velocity
Vol_M += - ufl.inner(dar, ufl.grad(del_p)) * dx # darcy

weak_form = Mom_B + Vol_M

#-------------------------------Solver----------------------------#

# Initialise non-linear problem
problem = NonlinearProblem(weak_form, up, list_bc)

# Initialise newton solver
solver = NewtonSolver(domain.comm, problem)

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

# Configure mumps
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

#-----------------------------Solution Settings-----------------------------#

# Create solution variables
Vu_sol, up_to_u_sol = V.sub(0).collapse() # dofs in u func subspace
u_sol = fem.Function(Vu_sol)

Vp_sol, up_to_p_sol = V.sub(1).collapse() # dofs in u func subspace
p_sol = fem.Function(Vp_sol)

V_vec = fem.VectorFunctionSpace(domain,("CG", 2)) #
vel_s_o = fem.Function(u_sp)

V_tensor = fem.FunctionSpace(domain,("CG", 2, (3,3))) # 2 is for 2nd order elements. (3,3) is size of stress tensor
stress_o = fem.Function(V_tensor)

# Force top surface space and function
force_FS = fem.VectorFunctionSpace(domain,("CG", 1))  #fem.VectorFunctionSpace(domain, "Real", 0)
df = ufl.TestFunction(force_FS)

# Expressions for output
w_FS_out_exp = fem.Expression(w_FS, u_sp.element.interpolation_points())
stress_out_exp = fem.Expression(P, u_sp.element.interpolation_points())
vel_s_out_exp = fem.Expression(vel_s, u_sp.element.interpolation_points())

# Rename for Paraview file
u_sol.name = "displacement"
p_sol.name = "pressure"
vel_s_o.name = "vel_s"
stress_o.name = "stress"

# Writing data to files
file1 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"disp.bp", u_sol, "bp4")
file2 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"press.bp", p_sol, "bp4")
file4 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"vel_s.bp", vel_s_o, "bp4")
file5 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"stress.bp", stress_o, "bp4") #

#---------------------------Solve problem-------------------------------------#
duration_solve = 0.0

for n in range(num_steps):

u_bc_2.value = axial_strain(t)
print('disp =',u_bc_2.value)

# Solve newton-steps
log.set_log_level(log.LogLevel.INFO)
duration_solve -= MPI.Wtime()
niter, converged = solver.solve(up)
duration_solve += MPI.Wtime()

PETSc.Sys.Print(
"Phys. Time {:.4f}, Calc. Time {:.4f}, Num. Iter. {}".format(
t, duration_solve, niter
)
)

u_n.x.array[:] = up.x.array[up_to_u]

u_sl = up.sub(0).collapse()
p_sl = up.sub(1).collapse()

u_sl.x.scatter_forward()
p_sl.x.scatter_forward()

u_sol.interpolate(u_sl)
p_sol.interpolate(p_sl)
vel_s_o.interpolate(vel_s_out_exp)
stress_o.interpolate(stress_out_exp) #

#---------stress on top surface-------#
traction = ufl.inner(df, ufl.dot(P,normal_vec))
Top_force = fem.assemble_scalar(fem.form(traction * ds(2)))
Top_stress = Top_force / (3.14 * pow(0.105,2))
print('traction =',Top_stress)

file1.write(t)
file2.write(t)
file3.write(t)
file4.write(t)
file5.write(t)
file6.write(t) #

#Update time
t += dt

file1.close()
file2.close()
file3.close()
file4.close()
file5.close()
file6.close() #


The msh file usedcan be found on this link: cylinder_21.msh - Google Drive

You send in a tuple to locate_dofs_topological when you have to map degrees of freedom from a function of a collapsed space to the degrees of freedom on the non-collapsed sub space.

When sending in a Constant, that is not associated with a function space, you do not need to create this map, thus one does what I showed in the code above.

The traction will just be a scalar (constant) value, so I do not see why you would save on the mesh?

Hi Dokken,

Thanks for the explanation. Yes you are right. My mistake. I was thinking about nodal reactions actually but I figured it out. Thank you for your time