Thank you very much for your help. With these changes, it seems like my code is running.
Sorry about the outside code, I changed it in the new version of the code that you can find below.
Now I’m trying to get the new shape of my beam. I saw in your example here: Hermite functionspace with unitinterval mesh - #3 by mscroggs that I need to separate the rotation and displacement solutions. What I didn’t really understand is how, in my 2D case, I can get the x and y displacements of the beam. Do I need to do something like : x_disp = disp * cos(rot) and y_disp = disp * sin(rot) ? I’m not sure I fully grasp what uh.x.array means.
from dolfinx import mesh, fem, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import basix
import numpy as np
import matplotlib.pyplot as plt
# from tools import EI_calc_static, menger_curvature
########## Initial parameters ##########
l_span = 4.
diameter = 31.05e-3
radius = diameter / 2
A = np.pi * (radius ** 2)
V = A * l_span
E = 62e8
EI_max = 2155
EI_min = 28.3
g = 10
rho = 270
chi_0 = 1.
########## Geometry initialization ##########
dim = 2
nb_cells = 200
beam_element = basix.ufl.element(basix.ElementFamily.Hermite, basix.CellType.interval, 3)
xn = np.linspace(0, l_span, nb_cells + 1)
yn = 0.005 * ((xn - l_span / 2) ** 2 - l_span)
x = np.vstack((xn, yn)).T
cs = np.zeros((nb_cells, 2), dtype=np.int64)
for i in range(nb_cells):
cs[i, 0] = i
cs[i, 1] = i + 1
cells = np.array(cs, dtype=np.int64)
shape = 'interval'
degree = 1
cell = ufl.Cell(shape)
domain = mesh.create_mesh(MPI.COMM_WORLD,
cells,
x,
ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree)))
########## Mesh and domain ##########
V = fem.FunctionSpace(domain, beam_element)
fdim = domain.topology.dim - 1
########## Pre-processing bretelle shape ##########
# # Curvature
# chi_new = menger_curvature(x.T[0], x.T[1])
# # Bending stiffness depending on curvature
# EI_new = EI_calc_static(domain, EI_max, EI_min, chi_0, chi_new)
EI_new = fem.Constant(domain, 2000.)
##########
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
########## Boundary conditions ##########
ubc = fem.Function(V)
with ubc.vector.localForm() as uloc:
uloc.set(0.)
# locate endpoints
startpt = mesh.locate_entities_boundary(domain, 0, lambda x: np.isclose(x[0], 0))
endpt = mesh.locate_entities_boundary(domain, 0, lambda x: np.isclose(x[0], l_span))
# locate DOFs from endpoints
startdof = fem.locate_dofs_topological(V, 0, startpt)
enddof = fem.locate_dofs_topological(V, 0, endpt)
# fix displacement of start point and rotation as well
fixed_disp1 = fem.dirichletbc(ubc, [startdof[0]])
fixed_rot1 = fem.dirichletbc(ubc, [startdof[1]])
fixed_disp2 = fem.dirichletbc(ubc, [enddof[0]])
fixed_rot2 = fem.dirichletbc(ubc, [enddof[1]])
bcs = [fixed_disp1, fixed_rot1, fixed_disp2, fixed_rot2]
########## Forces and tensions ##########
# distributed load value (due to weight)
q = fem.Constant(domain, -rho * A * g)
########## Solver ##########
def M(u):
return EI_new * ufl.div(ufl.grad(u))
a = ufl.inner(ufl.div(ufl.grad(u)), M(v)) * ufl.dx
L = q * v * ufl.dx
# Solver
problem = LinearProblem(a, L, bcs=bcs)
uh = problem.solve()
########## Post-processing ##########
# NOTE: The solution uh contains both the rotation and the displacement solutions
# The rotation and displacement 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)
new_coords = np.zeros((len(disp), 2))
for i in range(len(disp)):
new_coords[i][0] = xn[i] + disp[i] * np.cos(rot[i])
new_coords[i][1] = yn[i] + disp[i] * np.sin(rot[i])
# Plotting
figure, axis = plt.subplots(1, 1)
axis.plot(xn, yn, label='initial')
axis.plot(new_coords.T[0], new_coords.T[1], marker='.', linestyle='-', label='final')
axis.set_xlabel('x (m)')
axis.set_ylabel('y (m)')
axis.set_title("Displacement of FEM")
axis.legend()
plt.show()
Thank you very much for your quick reply!