Consistent reaction forces with mixed element

General
1

Hello everyone,

I’m trying to compute consistent reaction forces using this tutorial. However, my problem involves a coupled field with mixed finite elements.

I found one related post and was able to implement a similar approach, but unfortunately, the reaction forces I’m getting are incorrect. I’m having trouble identifying where my implementation is going wrong. If someone can help me out…

This was my attempt to set up such a scenario:

#definition
P = basix.ufl.element("CG", dmesh.basix_cell(), 1, shape=(dmesh.geometry.dim,))
S = basix.ufl.element("CG", dmesh.basix_cell(), 1)
R = basix.ufl.element("DG", dmesh.basix_cell(), 0)
# mixed element
mix = basix.ufl.mixed_element([P,S])
de
V = fem.functionspace(dmesh,mix)
W = fem.functionspace(dmesh,R)
du , dd = ufl.TestFunctions(V)

sol = fem.Function(V)
solold = fem.Function(V)

u, d = ufl.split(sol)
uold, dold = ufl.split(solold)
n = ufl.FacetNormal(dmesh)
(...)
bot_dofs = fem.locate_dofs_topological((V.sub(0),V0), 1, bot_facets)
top_dofs = fem.locate_dofs_topological((V.sub(0).sub(0),V01), 1, top_facets)
top_dofs_y = fem.locate_dofs_topological((V.sub(0).sub(1),V02), 1, top_facets)
crack_dofs = fem.locate_dofs_topological((V.sub(1),V1), 1, crack_facets)
(...)
mech = ufl.inner(sigma(u, dold),epsilon(du))*ufl.dx - ufl.dot(b,du)*ufl.dx - ufl.dot(t_bar,du)*ds_top
frac = G_cr*l*ufl.dot(ufl.grad(d),ufl.grad(dd)) *ufl.dx + ((G_cr/l) + 2*H(uold,u,Hold))*ufl.inner(d,dd) * ufl.dx - 2*dd*H(uold,u,Hold) *ufl.dx
F = mech + frac
(...)
# calculation of reaction force at the traction boundary
    v_reac = fem.Function(V)
    virtual_work_form = dolfinx.fem.form(ufl.action(F, v_reac))
    u_bc = fem.Function(V)
    u_bc.sub(0).sub(0).interpolate(one)
    toppedge = [dolfinx.fem.dirichletbc(u_bc, top_dofs,V.sub(0).sub(0))] 
    dolfinx.fem.set_bc(v_reac.vector, toppedge)
    print(f"Horizontal reaction Rx = {fem.assemble_scalar(virtual_work_form):.6f}",flush=True)
    
    u_bc.vector.set(0.0)
    v_reac.vector.set(0.0)
    u_bc.sub(0).sub(1).interpolate(one)
    toppedge = [dolfinx.fem.dirichletbc(u_bc, top_dofs_y,V.sub(0).sub(1))] 
    dolfinx.fem.set_bc(v_reac.vector, toppedge)
    print(f"Vertical reaction Ry = {fem.assemble_scalar(virtual_work_form):.6f}",flush=True)
    
    tractions = ufl.dot(sigma(u,d),n)
    forcex = fem.assemble_scalar(fem.form(tractions[0]*ds_top))    
    forcey = fem.assemble_scalar(fem.form(tractions[1]*ds_top))
    print(f"forcex = {forcex}",flush=True)
    print(f"forcey = {forcey}",flush=True)

i got this as result

Horizontal reaction Rx = -0.000000
Vertical reaction Ry = 0.000000
forcex = 3.9484366369161937e-16
forcey = 4.947636518989789e-14

Horizontal reaction Rx = 31578.132069
Vertical reaction Ry = 8781.043925
forcex = 92977.05888077534
forcey = -732.5915813075735

Horizontal reaction Rx = 63146.070821
Vertical reaction Ry = 17415.128594
forcex = 185818.18618585676
forcey = -1843.9803530233448