I want to solve an elastic cantilever beam under compression problem.The left side of the beam is fixed and the upper end is subjected to positive stresses, considering the effect of gravity.The code is as follows：
from fenics import *
import numpy as np
def solver_bcs(boundary_conditions,
degree=1,
subdomains=[],
linear_solver='Krylov',
abs_tol=1E-5,
rel_tol=1E-3,
max_iter=1000):
# Scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma
# Define strain and stress
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
# mesh = UnitSquareMesh(Nx, Ny)
V = FunctionSpace(mesh, 'P', degree)
# Define boundary subdomains
tol = 1e-14
class BoundaryX0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol)
class BoundaryX1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 10, tol)
class BoundaryY0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0, tol)
class BoundaryY1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 3, tol)
class BoundaryZ0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0, tol)
class BoundaryZ1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 3, tol)
# Mark boundaries
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
boundary_markers.set_all(9999)
bx0 = BoundaryX0()
bx1 = BoundaryX1()
by0 = BoundaryY0()
by1 = BoundaryY1()
bz0 = BoundaryZ0()
bz1 = BoundaryZ1()
bx0.mark(boundary_markers, 0)
bx1.mark(boundary_markers, 1)
by0.mark(boundary_markers, 2)
by1.mark(boundary_markers, 3)
bz0.mark(boundary_markers, 4)
bz1.mark(boundary_markers, 5)
# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
boundary_markers, i)
bcs.append(bc)
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
# Collect Neumann integrals
integrals_N = []
for i in boundary_conditions:
if 'Neumann' in boundary_conditions[i]:
if boundary_conditions[i]['Neumann'] != 0:
g = boundary_conditions[i]['Neumann']
integrals_N.append(dot(T, v)*ds(i))
# Sum integrals to define variational problem
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + sum(integrals_N)
# Set linear solver parameters
prm = LinearVariationalSolver.default_parameters()
if linear_solver == 'Krylov':
prm["linear_solver"] = 'gmres'
prm["preconditioner"] = 'ilu'
prm["krylov_solver"]["absolute_tolerance"] = abs_tol
prm["krylov_solver"]["relative_tolerance"] = rel_tol
prm["krylov_solver"]["maximum_iterations"] = max_iter
else:
prm["linear_solver"] = 'lu'
# Compute solution
u = Function(V)
solve(a == L, u, bcs, solver_parameters=prm)
return u
def demo_bcs():
"Compute and plot solution using a combination of boundary conditions"
# Define manufactured solution in sympy and derive f, g, etc.
import sympy as sym
x, y = sym.symbols('x[0], x[1]') # needed by UFL
u_0 = Constant((0, 0, 0))
T_1 = Constant((0, 0, 0))
T_2 = Constant((0, 0, 0))
T_3 = Constant((0, 1, 0))
T_4 = Constant((0, 0, 0))
T_5 = Constant((0, 0, 0))
# # # Collect variables
# # variables = [u_0, T_1, T_2, T_3, T_4, T_5]
# # # Turn into C/C++ code strings
# # variables = [sym.printing.ccode(var) for var in variables]
# # # Turn into FEniCS Expressions
# # variables = [Expression(var, degree=3) for var in variables]
# # Extract variables
# u_0, T_1, T_2, T_3, T_4, T_5 = variables
# Define boundary conditions
boundary_conditions = {0: {'Dirichlet': u_0}, # x = 0
1: {'Neumann': T_1}, # x = 1
2: {'Neumann': T_2}, # y = 0
3: {'Neumann': T_3}, # y = 3
4: {'Neumann': T_4}, # z = 0
5: {'Neumann': T_5}} # z = 3
# Compute solution
u = solver_bcs(boundary_conditions,
degree=1, linear_solver='direct')
# Save and plot solution
vtkfile = File('/solution_bcs.pvd')
vtkfile << u
# plot(u)
demo_bcs()
However, the following error is reported at runtime:
RuntimeError Traceback (most recent call last)
Cell In[4], line 170
167 vtkfile << u
168 # plot(u)
--> 170 demo_bcs()
Cell In[4], line 162, in demo_bcs()
154 boundary_conditions = {0: {'Dirichlet': u_0}, # x = 0
155 1: {'Neumann': T_1}, # x = 1
156 2: {'Neumann': T_2}, # y = 0
157 3: {'Neumann': T_3}, # y = 3
158 4: {'Neumann': T_4}, # z = 0
159 5: {'Neumann': T_5}} # z = 3
161 # Compute solution
--> 162 u = solver_bcs(boundary_conditions,
163 degree=1, linear_solver='direct')
165 # Save and plot solution
166 vtkfile = File('/solution_bcs.pvd')
Cell In[4], line 85, in solver_bcs(boundary_conditions, degree, subdomains, linear_solver, abs_tol, rel_tol, max_iter)
83 for i in boundary_conditions:
84 if 'Dirichlet' in boundary_conditions[i]:
---> 85 bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'],
86 boundary_markers, i)
87 bcs.append(bc)
89 # Define trial and test functions
File ~/anaconda3/envs/fenicsproject/lib/python3.12/site-packages/dolfin/fem/dirichletbc.py:131, in DirichletBC.__init__(self, *args, **kwargs)
128 if (len(kwargs) > 0):
129 raise RuntimeError("Invalid keyword arguments", kwargs)
--> 131 super().__init__(*args)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to create Dirichlet boundary condition.
*** Reason: Expecting a scalar boundary value but given function is vector-valued.
*** Where: This error was encountered inside DirichletBC.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 2e001bd1aae8e14d758264f77382245e6eed04b0
*** -------------------------------------------------------------------------