Getting a recursion error for DirichletBC

General
1

Hi eveyone, I am new to fenics, I am trying to solve a cantilever beam problem with point load at free end. However I am getting following error when I run the code.

Traceback (most recent call last):
File “new.py”, line 92, in
main()
File “new.py”, line 36, in main
bc = DirichletBC(V, Constant((0, 0)), root_boundary)
File “/usr/local/lib/python3.6/dist-packages/ufl/coefficient.py”, line 125, in Constant
domain = as_domain(domain)
File “/usr/local/lib/python3.6/dist-packages/ufl/domain.py”, line 296, in as_domain
cell = as_cell(domain)
File “/usr/local/lib/python3.6/dist-packages/ufl/cell.py”, line 327, in as_cell
return TensorProductCell(cell)
File “/usr/local/lib/python3.6/dist-packages/ufl/cell.py”, line 223, in init
self._cells = tuple(as_cell(cell) for cell in cells)
File “/usr/local/lib/python3.6/dist-packages/ufl/cell.py”, line 223, in
self._cells = tuple(as_cell(cell) for cell in cells)
File “/usr/local/lib/python3.6/dist-packages/ufl/cell.py”, line 327, in as_cell
return TensorProductCell(cell)

RecursionError: maximum recursion depth exceeded while calling a Python object

here is my code:

from dolfin import *
from ufl import *
import matplotlib.pyplot as plt


# Cantilever beam with point load

L = 0.5 # length
b = 0.05 # width
E = 70*(10**9) # modulus of elasticity
nu = 0.3 # poissons ratio
Lambda = (E*nu)/((1+nu)*(1-2*nu)) # lames parameter
Mu = E/(2*(1+nu)) # lames parameter

def main():

	# Meshing and creating function space	

	nx = 50 # no. of cells along length
	ny = 5 # no. of cells along width

	Mesh = RectangleMesh(Point(0.0, 0.0), Point(L, b), nx, ny)
	plot(Mesh)
	plt.savefig('Mesh.png')

	V = VectorFunctionSpace(Mesh, 'Lagrange', 1)


	# Boundary conditons
	def root_boundary(x, on_boundary):
		tol = 1e-14
		return on_boundary and x[0]<tol
		
	bc = DirichletBC(V, Constant((0, 0)), root_boundary)
	
	# Application of point load
	
	class PointLoad(UserExpression):

		def __init__(self, **kwargs):
		    super().__init__(degree=kwargs["degree"])

		    self.point = kwargs['point']
		    self.value = kwargs['value']

		def eval(self, value, x):

		    if near(x[0], self.point[0]) and near(x[1], self.point[1]):
		        value[0] = self.value[0]
		        value[1] = self.value[1]
		    else:  
		        value[0] = 0
		        value[1] = 0

		def value_shape(self):
		    return (2,)
		    

	# Stress and Strain Tensors

	def small_strain(u):
		return sym(nabla_grad(u))
	
	def cauchy_stress(u):
		return Lambda*(tr(small_strain(u)))*Identity(2) + 2*Mu*small_strain(u)
	
	# Variational form
		
	u_trail = TrailFunction(V)
	v_test = TestFunction(V)
	f = PointLoad(point=(L, b), value=(0,-10000), degree=1)
	
	lhs = inner(cauchy_stress(u_trail), small_strain(v_test))*dx
	rhs = dot(f, v_test)*dx
	
	# Solving
	u_solution = Function(V)
	solve(lhs==rhs, u_solution, bc)
	plot(u_solution)
	plt.savefig('u_field_image.png')
	u_magnitude = sqrt(dot(u, u))
	u_magnitude = project(u_magnitude, V)
	plot(u_magnitude, 'Displacement magnitude')
	plt.savefig('u_magnitude.png')
	print('min/max u:',
		  u_magnitude.vector().get_local().min(),
		  u_magnitude.vector().get_local().max())

if __name__ == "__main__":
    main()

Thank you in advance

2

Please consider the following patch to your file

diff --git a/new.py b/new.py
index e7ada69..151927f 100644
--- a/new.py
+++ b/new.py
@@ -1,5 +1,5 @@
-from dolfin import *
 from ufl import *
+from dolfin import *
 import matplotlib.pyplot as plt
 
 
@@ -66,7 +66,7 @@ def main():
 
        # Variational form
 
-       u_trail = TrailFunction(V)
+       u_trail = TrialFunction(V)
        v_test = TestFunction(V)
        f = PointLoad(point=(L, b), value=(0,-10000), degree=1)
 
@@ -78,8 +78,9 @@ def main():
        solve(lhs==rhs, u_solution, bc)
        plot(u_solution)
        plt.savefig('u_field_image.png')
-       u_magnitude = sqrt(dot(u, u))
-       u_magnitude = project(u_magnitude, V)
+       u_magnitude = sqrt(dot(u_solution, u_solution))
+       scalar_V = FunctionSpace(Mesh, 'Lagrange', 1)
+       u_magnitude = project(u_magnitude, scalar_V)
        plot(u_magnitude, 'Displacement magnitude')
        plt.savefig('u_magnitude.png')
        print('min/max u:',

Apart from a couple of easily fixable typos, the key change is to swap the order of the imports (ufl first, than dolfin).

3

I would in general suggest never to use wildcard imports, due to the issues you see here with namespace clashes. This is also mentioned in: 6. Modules — Python 3.12.0 documentation

Although certain modules are designed to export only names that follow certain patterns when you use import *, it is still considered bad practice in production code.

4

Thank you francesco for the response.
Now my code is running without excution error.
However the displacement at the free end after running is not matching the analytical solution
I have considered a beam of following dimensions
length = 1m
width = 0.1m
Elasticity modulus (E) = 70*(10^9) N/m^2
poissons ratio (nu) = 0.3
area moment of inertia (I) = (width^4)/12
load at end (P) = -10000N
I have used the following formula for calculating end deflection
displacement(at free end) = (PL**3)/(3E*I)

I guess there is something wrong with lines(in code) concerning with application of point load at free end. If it was so, could you please suggest me a correct way of doing it. Feel happy to receive any other suggestions related to this problem. Thank you

here is the modified code

from ufl import *
from dolfin import *
import matplotlib.pyplot as plt


# Cantilever beam with point load

L = 1.0 # length
b = 0.1 # width
Mu = 26.92*(10**9)
Lambda = 40.384*(10**9)


def main():

	# Meshing and creating function space	

	nx = 50 # no. of cells along length
	ny = 5 # no. of cells along width

	Mesh = RectangleMesh(Point(0.0, 0.0), Point(L, b), nx, ny)
	plot(Mesh)
	plt.savefig('Mesh.png')

	V = VectorFunctionSpace(Mesh, 'Lagrange', 1)


	# Boundary conditons
	def root_boundary(x, on_boundary):
		tol = 1e-14
		return on_boundary and x[0]<tol
		
	bc = DirichletBC(V, Constant((0, 0)), root_boundary)
	
	# Application of point load
	
	class PointLoad(UserExpression):

		def __init__(self, **kwargs):
		    super().__init__(degree=kwargs["degree"])

		    self.point = kwargs['point']
		    self.value = kwargs['value']

		def eval(self, value, x):

		    if near(x[0], self.point[0]) and near(x[1], self.point[1]):
		        value[0] = self.value[0]
		        value[1] = self.value[1]
		    else:  
		        value[0] = 0
		        value[1] = 0

		def value_shape(self):
		    return (2,)
		    

	# Stress and Strain Tensors

	def small_strain(u):
		return sym(nabla_grad(u))
	
	def cauchy_stress(u):
		return Lambda*(tr(small_strain(u)))*Identity(2) + 2*Mu*small_strain(u)
	
	# Variational form
		
	u_trail = TrialFunction(V)
	v_test = TestFunction(V)
	f = PointLoad(point=(L, b), value=(0,-10000), degree=1)
	
	lhs = inner(cauchy_stress(u_trail), small_strain(v_test))*dx
	rhs = dot(f, v_test)*dx
	
	# Solving
	u_solution = Function(V)
	solve(lhs==rhs, u_solution, bc)
	print(abs(u_solution.vector().get_local().min()))
	plot(u_solution)
	plt.savefig('u_field_image.png')
	u_magnitude = sqrt(dot(u_solution, u_solution))
	scalar_V = FunctionSpace(Mesh, 'Lagrange', 1)
	u_magnitude = project(u_magnitude, scalar_V)
	plot(u_magnitude, 'Displacement magnitude')
	plt.savefig('u_magnitude.png')
	print((10000*12)/(3*70*(10**9)*(0.1**4)))

if __name__ == "__main__":
    main()
5

The behavior of the solution from the plots you save to file looks good to me, in the sense that the solution is zero on the left, varies linearly in the domain, and has maximum on the right-end side. The vectors are pointing downwards as expected.

The mesh is of course quite coarse, but that does not justify the values being different by orders of magnitude.

I would double check the formula you use for computing the exact solution, especially if all the terms appearing there are dimensionally consistent (e.g., some lengths are measured in m while others in cm, or similar).

For the comparison with the numerical solution, you should use

print(u_magnitude.vector().get_local().min(), u_magnitude.vector().get_local().max())

The print you have

print(abs(u_solution.vector().get_local().min()))

is taking the absolute value of a component (either x or y) of the displacement, rather than the norm of the displacement vector.