Hi everyone,
I’m working on a 2D elasticity problem in FEniCS where I need to apply a Neumann boundary condition to the boundary of a circular cavity within my mesh, i’m using fenics 2019.1.0 version.
My Setup:
Mesh: A rectangular mesh with a circular cavity at its center.
Function Space: ‘VectorFunctionSpace’ for displacement.
Subdomain & Boundary Marking:
Cells within the cavity are marked using ‘MeshFunction’.
The boundary of this cavity is identified and marked to apply a Neumann condition.
**Boundary Condition:**I’m trying to apply a traction force ('T_cavity = 100 * n`) on the cavity boundary using a custom ‘ds’ measure.
Problem:
After running the simulation, the displacement field returns zero everywhere. This suggests that the Neumann condition isn’t being applied correctly, even though the boundary markers seem fine. here is my code wich recreates the issue:
from dolfin import *
from fenics import *
import numpy as np
#====mecanical parameters
nx, ny = 100, 100
E = 100000
nu = 0.3
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
#==== mesh and vector function space definition
mesh = RectangleMesh(Point(0, 0), Point(nx, ny), nx, ny, 'left/right')
Vvector = VectorFunctionSpace(mesh, "P", 2)
#==== useful function
def eps(u):
return sym(grad(u))
def sigma(u):
return lmbda*div(u)*Identity(2) + 2*mu*eps(u)
#==== defining dirichlet boundary condition on the domain boundary
def add_boundary_conditions(FunctionSpace, boundary_condition):
class Boundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
boundary = Boundary()
bc = DirichletBC(FunctionSpace, Constant(boundary_condition), boundary)
return [bc]
#====this part is the code snippet used to reproduce the conditions which represent how mycavity subdomain is defined
#==function to mark cell which belongs to the cavity
def mark_subdomain_cells():
subdomain_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
circle_center = Point(50, 50) # Define the center of the circle
radius = 25 # Define the radius of the circle
for cell in cells(mesh):
cell_center = cell.midpoint()
if cell_center.distance(circle_center) <= radius:
subdomain_markers[cell] = 1
return subdomain_markers
#==funtion to mark with tag 1 the exterior facets of the cavity to define the cavity boundary
def mark_boundary_facets(subdomain_markers):
# Initialize connectivity needed to find neighboring cells
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1, 0)
mesh.init(mesh.topology().dim() - 1, mesh.topology().dim())
for cell in cells(mesh):
if subdomain_markers[cell.index()] == 1:
for facet in facets(cell):
is_boundary_facet = False
facet_cells = facet.entities(mesh.topology().dim())
if len(facet_cells) == 1:
# Facet is on the boundary of the mesh
is_boundary_facet = True
elif len(facet_cells) == 2:
# Facet is internal, check if it separates cells with different markers
neighbor_cell_index = facet_cells[0] if facet_cells[1] == cell.index() else facet_cells[1]
if subdomain_markers[neighbor_cell_index] != 1:
is_boundary_facet = True
if is_boundary_facet:
boundary_markers[facet.index()] = 1
return boundary_markers
#==run simulation
def run_finite_elements():
# Define the variational problem
u = TrialFunction(Vvector)
v = TestFunction(Vvector)
n = FacetNormal(mesh)
T_cavity = 100 * n
boundary_conditions = add_boundary_conditions(Vvector,(0,0))
subdomain_markers = mark_subdomain_cells()
boundary_markers = mark_boundary_facets(subdomain_markers)
ds1 = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
a = inner(sigma(u), eps(v)) * dx
L = dot(T_cavity, v) * ds1(1)
u = Function(Vvector)
solve(a == L, u, boundary_conditions)
print("Max displacement =", np.max(u.vector().get_local()))
print("Min displacement =", np.min(u.vector().get_local()))
File("boundary.pvd") << boundary_markers
File("subdomain.pvd") << subdomain_markers
File("displacement.pvd") << u
run_finite_elements()