How to define a spring boundary condition on euler bernoulli beam theory?

Hi,
I’m trying to set my boundary conditions as a spring in a 2D space with 1D wire using euler - Bernoulli beam theory. Can someone help me with this probably very easy question? I’m just started learning fenicsx and I’m already stuck.

I have tried the following code which I found in the forum and tried to modify to my code:

``````tag = meshtags(domain, 1, startpt, [1])
ds = Measure("ds", domain=domain, subdomain_data=tag)
k_form += k_spring * u_ * v * ds(0)
``````

But this code return an error which I don’t know how to solve:

``````AssertionError
Exception ignored in: <function LinearProblem.__del__ at 0x7fd9272668e0>
Traceback (most recent call last):
File "/home/maarten/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 805, in __del__
self._solver.destroy()
^^^^^^^^^^^^
AttributeError: 'LinearProblem' object has no attribute '_solver'
``````

Here is my working code with dirichlet BC:

``````import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, dirichletbc,functionspace,
locate_dofs_topological,Constant)
from dolfinx.io import VTKFile
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities, locate_entities_boundary,create_interval
from ufl import (SpatialCoordinate,inner, TestFunction, TrialFunction, div, grad, dx)

E = 210e9 # N/m^2
L = 10.0 # m
thickness = 0.005 # m
h = 0.12 # m
b = 0.09 # m
rho= 7.8e-3
g = 9.81
k_spring = 4000 / 0.004 # N/m
f = 1000.0 # N/m
NUM_ELEM = int(L * 100)
domain = create_interval(MPI.COMM_WORLD, NUM_ELEM, [0, L])
x = SpatialCoordinate(domain)
E = Constant(domain,E)
rho = Constant(domain,rho)
g = Constant(domain,g)
f = Constant(domain,f)
A = b*h
I = b*h**3/12
EI = E * I
EA = E * A
number_of_supports = round(L / 10)
supports = np.round(np.linspace(0, L, number_of_supports + 1), 2)

def M(u):

beam_element=basix.ufl.element(basix.ElementFamily.Hermite, basix.CellType.interval, 3)

W = functionspace(domain,beam_element)
u_ = TestFunction(W)
v = TrialFunction(W)

#bilinear form (LHS)
#linear form construction (RHS)
l_form = f*u_*dx

ubc = Function(W)
with ubc.vector.localForm() as uloc:
uloc.set(0.)

# list of boundary conditions
bcs = []
for pt in supports:
startpt=locate_entities(domain,0,lambda x : np.isclose(x[0], pt))
#locate DOFs from pt
dof=locate_dofs_topological(W,0,startpt)
#fix displacement of start point and rotation as well
bcs.append(dirichletbc(ubc,np.array([dof[0]])))

u = Function(W)

problem = LinearProblem(k_form, l_form, u=u, bcs=bcs)
uh=problem.solve()
``````

The actually relevant error probably is

``````Traceback (most recent call last):
File "/tmp/a.py", line 58, in <module>
dof=locate_dofs_topological(W,0,startpt)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "lib/python3.11/site-packages/dolfinx/fem/bcs.py", line 87, in locate_dofs_topological
return _cpp.fem.locate_dofs_topological(V._cpp_object, entity_dim, _entities, remote)  # type: ignore
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Entity-to-cell connectivity has not been computed.
``````

That error you get can be fixed by adding a `domain.topology.create_connectivity(0, domain.topology.dim)` before the final `for` loop.

1 Like

I have added the line but now it throws another error.

``````import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, dirichletbc,functionspace,
locate_dofs_topological,Constant)
from dolfinx.io import VTKFile
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import locate_entities, locate_entities_boundary,create_interval, meshtags
from ufl import (SpatialCoordinate,inner, TestFunction, TrialFunction, div, grad, dx, Measure)

E = 210e9 # N/m^2
L = 10.0 # m
thickness = 0.005 # m
h = 0.12 # m
b = 0.09 # m
rho= 7.8e-3
g = 9.81
k_spring = 4000 / 0.004 # N/m
f = 1000.0 # N/m
NUM_ELEM = int(L * 100)
domain = create_interval(MPI.COMM_WORLD, NUM_ELEM, [0, L])
x = SpatialCoordinate(domain)
E = Constant(domain,E)
rho = Constant(domain,rho)
g = Constant(domain,g)
f = Constant(domain,f)
A = b*h
I = b*h**3/12
EI = E * I
EA = E * A
number_of_supports = round(L / 10)
supports = np.round(np.linspace(0, L, number_of_supports + 1), 2)

def M(u):

beam_element=basix.ufl.element(basix.ElementFamily.Hermite, basix.CellType.interval, 3)

W = functionspace(domain,beam_element)
u_ = TestFunction(W)
v = TrialFunction(W)

#bilinear form (LHS)
#linear form construction (RHS)
l_form = f*u_*dx

ubc = Function(W)
with ubc.vector.localForm() as uloc:
uloc.set(0.)

# list of boundary conditions
bcs = []
domain.topology.create_connectivity(0, domain.topology.dim)
for pt in supports:
startpt=locate_entities(domain,0,lambda x : np.isclose(x[0], pt))
#locate DOFs from pt
dof=locate_dofs_topological(W,0,startpt)
#fix displacement of start point and rotation as well
tag = meshtags(domain, 1, startpt, [1])
ds = Measure("ds", domain=domain, subdomain_data=tag)
k_form += k_spring * u_ * v * ds(0)

u = Function(W)

problem = LinearProblem(k_form, l_form, u=u, bcs=bcs)
uh=problem.solve()
``````

error:

``````AssertionError
Exception ignored in: <function LinearProblem.__del__ at 0x7f4b3b4c2840>
Traceback (most recent call last):
File "/home/maarten/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 805, in __del__
self._solver.destroy()
^^^^^^^^^^^^
AttributeError: 'LinearProblem' object has no attribute '_solver'
``````

Can’t reproduce. You’ll need to tell us exactly how you installed `dolfinx`, and which version you are running.

Note however that

so it’s probably safe to ignore that exception.

Did you try the code from my second post? The code in my original post uses dirichletBC instead of springs and that one works on my computer as well.
I installed fenics using pip in a conda environment.
Python 3.12.2
fenics-basix 0.8.0
fenics-dolfinx 0.8.0
fenics-ffcx 0.8.0
fenics-ufl 2024.1.0

I modified the code a bit. Now it does run and the spring is correctly setup at the first and last node but not in nodes in the middle. What am I doing wrong?

Part of the code I modified:

``````domain.topology.create_connectivity(0, domain.topology.dim)
pts = np.array([], dtype=np.int32)
for pt in supports:
startpt=locate_entities(domain,0,lambda x : np.isclose(x[0], pt))
pts = np.concatenate((pts, startpt))

tag = meshtags(domain, int(domain.topology.dim - 1), pts, np.ones_like(pts))
ds = Measure("ds", domain=domain, subdomain_data=tag)
k_form += k_spring * inner(v, u_) * ds(1)
``````

Full code:

``````import numpy as np
from mpi4py import MPI
import basix
from dolfinx.fem import (Function, dirichletbc,functionspace, form, set_bc, assemble_scalar,
locate_dofs_topological,Constant,Expression)
from dolfinx.io import VTKFile
from dolfinx.fem.petsc import LinearProblem, assemble_matrix, assemble_vector, apply_lifting
from dolfinx.mesh import locate_entities, locate_entities_boundary,create_interval, meshtags
from ufl import (SpatialCoordinate,inner, dot, TestFunction, TrialFunction, action, div, grad, dx, Measure, derivative, conditional, le, lhs, rhs)
import matplotlib.pyplot as plt
from section import inertia, area
import dolfinx.fem.petsc as petsc
import time

t1 = time.time()

def last_value(arr):
if len(arr) == 0:
return -1
return arr[-1]

#################################################################
########### SUPPORTS PARAMETERS #################################
#################################################################
supports = [0.0,5.0,10.0]
k_spring = 4000 / 0.004

########## GEOMETRIC INPUT ####################
E = 210e9 # N/m^2
L = 10.0 # m
thickness = 0.005 # m
h = 0.12 # m
b = 0.09 # m

rho= 7.8e-3
g = 9.81

#################################################################
########### CONSTRUCT BEAM MESH #################################
#################################################################
NUM_ELEM = int(L * 100)
domain = create_interval(MPI.COMM_WORLD, NUM_ELEM, [0, L])

#################################################################
##### ENTER MATERIAL PARAMETERS AND CONSTITUTIVE MODEL ##########
#################################################################
x = SpatialCoordinate(domain)

E = Constant(domain,E)
rho = Constant(domain,rho)
g = Constant(domain,g)
k_spring = Constant(domain,k_spring)
f = Constant(domain,1000.0)

I = b * h ** 3 / 12
EI = E * I

#################################################################
########### COMPUTE STATIC SOLUTION #############################
#################################################################

#define Moment expression
def M(u):

# Create Hermite order 3 on a interval (for more informations see:
#    https://defelement.com/elements/examples/interval-Hermite-3.html )
beam_element=basix.ufl.element(basix.ElementFamily.Hermite, basix.CellType.interval, 3)

#finite element function space on domain, with trial and test fxns
W = functionspace(domain,beam_element)
print("Number of DOFs: %d" % W.dofmap.index_map.size_global)
print("Number of elements (intervals): %d" % NUM_ELEM)
print("Number of nodes: %d" % (NUM_ELEM+1))
u_ = TestFunction(W)
v = TrialFunction(W)

#bilinear form (LHS)

#linear form construction (RHS)
l_form = f*u_*dx

#APPLY BOUNDARY CONDITIONS
#initialize function for boundary condition application
ubc = Function(W)
with ubc.vector.localForm() as uloc:
uloc.set(0.)

# list of boundary conditions
bcs = []
domain.topology.create_connectivity(0, domain.topology.dim)
pts = np.array([], dtype=np.int32)
for pt in supports:
startpt=locate_entities(domain,0,lambda x : np.isclose(x[0], pt))
pts = np.concatenate((pts, startpt))

tag = meshtags(domain, int(domain.topology.dim - 1), pts, np.ones_like(pts))
ds = Measure("ds", domain=domain, subdomain_data=tag)
k_form += k_spring * inner(v, u_) * ds(1)

#SOLVE VARIATIONAL PROBLEM
#initialize function in functionspace for beam properties
u = Function(W)

# solve variational problem
problem = LinearProblem(k_form, l_form, u=u, bcs=bcs)
uh=problem.solve()
#uh.name = "Displacement and Rotation"

# Define the shear force expression
shear_force_expr = EI * uh_xxx

V = functionspace(domain, ("DG", 0))
shear_expr = Expression(shear_force_expr, V.element.interpolation_points())
shear = Function(V)
shear.interpolate(shear_expr)

reaction = []
for support in supports:
index=locate_entities(domain,0,lambda x : np.isclose(x[0], support))[0]

if index == 0:
left_side = 0
else:
left_side = shear.x.array[index-1]

if index == len(shear.x.array):
right_side = 0
else:
right_side = shear.x.array[index]
f_reac = right_side - left_side
reaction.append(f_reac)
print(f"Support at {support} has reaction force of {round(f_reac / 1000, 2)} kN")

print("time taken: ", time.time()-t1)
#################################################################
########### SAVE AND VISUALIZE RESULTS ##########################
#################################################################

#NOTE: The solution uh contains both the rotation and the displacement solutions
#The rotation and displacment solutions can be separated as follows:
disp = np.empty(0)
rot = np.empty(0)
for i,x in enumerate(uh.x.array):
if i % 2 != 0:
rot = np.append(rot,x)
else:
disp = np.append(disp,x)

x_fem = np.linspace(0,L,disp.shape[0])

figure, axis = plt.subplots(1,1)

print("Maximum magnitude displacement is: %e" % np.min(disp))
``````