[BUG]: Non-interacting cylinders bend in unintuitive direction.

Summarize the issue

I have two cylinders which are not supposed to interact with each other. I apply a pressure at the right face of both cylinders. But instead of going straight they are bending in a random direction.

If I move them further away from each other this effect increases. cylinder_extrem Moving the cylinders to different coordinates does not change the direction of the bending, e.g. if their positions were swapped they would bend towards each other and not away from each other.

I am not sure if this is a numerical artifact or if there is something not working as expected when there are multiple objects.

It is possible to reduce this effect by lowering the mesh size with gmsh.option.setNumber("Mesh.MeshSizeMax", 0.05), but I would expect it to also work with a bigger mesh as long as it works when there is only one cylinder (with the same parameter). single_cylinder

How to reproduce the bug

Run the provided code and look at the generated "results.pvd" using paraview. Changing plot_both_cylinders = False will only generate one cylinder.

Minimal Example (Python)

import numpy as np  # type: ignore
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction  # type: ignore
from mpi4py import MPI  # type: ignore
from dolfinx import fem, io, mesh
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.fem import locate_dofs_geometrical, DirichletBC
from petsc4py import PETSc # type: ignore
import gmsh

plot_both_cylinders = True

pressure_value = 8.0
E = 10 # Young's modulus
nu = 0 # Poisson's ratio
radius = 1
radius_2 = 1
length = 10

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gdim = 3

# add first cylinder
gmsh.model.add("cylinder")
cylinder_tag = gmsh.model.occ.addCylinder(0, 0, 0, length, 0, 0, radius)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(gdim, [cylinder_tag], 0)

surfaces = gmsh.model.getBoundary([(3, cylinder_tag)], oriented=False, recursive=False)
right_face_tag = surfaces[1][1]
gmsh.model.addPhysicalGroup(gdim - 1, [right_face_tag], 1)
gmsh.model.setPhysicalName(gdim - 1, 1, "Right")

left_face_tag = surfaces[0][1]
gmsh.model.addPhysicalGroup(gdim - 1, [left_face_tag], 2)
gmsh.model.setPhysicalName(gdim - 1, 2, "Left")

lateral_surface_tags = [surf[1] for surf in surfaces[2:]]
gmsh.model.addPhysicalGroup(gdim - 1, lateral_surface_tags, 3)
gmsh.model.setPhysicalName(gdim - 1, 3, "Lateral")

if plot_both_cylinders:
    # add second cylinder
    second_cylinder_tag = gmsh.model.occ.addCylinder(0, 0,20*radius,  length, 0, 0, radius_2)
    gmsh.model.occ.synchronize()
    gmsh.model.addPhysicalGroup(gdim, [second_cylinder_tag], 1)

    surfaces = gmsh.model.getBoundary([(3, second_cylinder_tag)], oriented=False, recursive=False)
    second_right_face_tag = surfaces[1][1]
    gmsh.model.addPhysicalGroup(gdim - 1, [second_right_face_tag], 4)
    gmsh.model.setPhysicalName(gdim - 1, 4, "secondRight")

    second_left_face_tag = surfaces[0][1]
    gmsh.model.addPhysicalGroup(gdim - 1, [second_left_face_tag], 5)
    gmsh.model.setPhysicalName(gdim - 1, 5, "secondLeft")

    second_lateral_surface_tags = [surf[1] for surf in surfaces[2:]]
    gmsh.model.addPhysicalGroup(gdim - 1, second_lateral_surface_tags, 6)
    gmsh.model.setPhysicalName(gdim - 1, 6, "secondLateral")

gmsh.model.mesh.generate(gdim)
domain, _, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=gdim)
gmsh.finalize()

shape = (gdim,)

# Define the function space
V = fem.functionspace(domain, ("P", 2, shape))
u = TrialFunction(V)
v = TestFunction(V)

u_h = fem.Function(V, name="Displacement")

lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

# Define strain and stress tensors
def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return lambda_ * tr(epsilon(u)) * Identity(len(u)) + 2 * mu * epsilon(u)

pressure = fem.Constant(domain, np.array([pressure_value, 0, 0]))  

dx = Measure("dx", domain=domain)       # Volume measure
ds = Measure("ds", domain=domain, subdomain_data=facets)  # Surface measure for boundaries

right_face = 1
right_face_second = 4

# Define the boundary conditions
def left(x):
    return np.isclose(x[0], 0)

left_dofs = locate_dofs_geometrical(V, left)

bcs = [
    fem.dirichletbc(np.zeros(3,), left_dofs, V) # keep the left face fixed
]

# Define the weak form
a = (inner(sigma(u), epsilon(v)) * dx)

if plot_both_cylinders:
    L = inner(pressure, v) * (ds(right_face) + ds(right_face_second))
else:
    L = inner(pressure, v) * ds(right_face)

bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form, bcs=bcs)
A.assemble()
b = create_vector(linear_form)

solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

problem = fem.petsc.LinearProblem(
    a, L, u=u_h, bcs=bcs,
    petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)

vtk = io.VTKFile(domain.comm, "results.pvd", "w")

# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
    loc_b.set(0)
assemble_vector(b, linear_form)

# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)

# Solve linear problem
solver.solve(b, u_h.x.petsc_vec)
u_h.x.scatter_forward()
vtk.write_function(u_h)

vtk.close()

Output (Python)

No response

Version

main branch

DOLFINx git commit

4116cca8cf91f1a7e877d38039519349b3e08620

Installation

Docker image

Additional information

No response

FEniCS / dolfinx