Create mesh for Euler-Bernoulli beam with Hermite elements

mesh
1

Hi everyone,
I’m trying to model a 1D wire in 2D space with fenicsx using Euler Bernoulli. Thus I need to have Hermite elements for my wire. I managed to create these elements using basix. Now I am having trouble creating a mesh for the wire with this type of element.

Here’s my code:

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
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)

# Option 1?
domain_init = ufl.Mesh(beam_element)
# Option 2?
# gdim, shape, degree = dim, "interval", 3
# cell = ufl.Cell(shape, geometric_dimension=gdim)
# domain_init = ufl.Mesh(ufl.VectorElement("CG", cell, degree))

xn = np.linspace(0, l_span, nb_cells+1)
yn = 0.5*((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)

########## Mesh and domain ##########

domain = mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain_init)
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)

##########

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))

# 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()

Any help would be greatly appreciated, thank you very much!

2

The mesh creation does not depend on the creation of the ufl element. The creation of a FunctionSpace is the first place that the finite element type (e.g. Hermite elements) and the mesh are connected.

In the case of a topologically 1D and geometrically 1D beam, you should be able to use the normal dolfinx syntax to create an interval mesh.

domain = dolfinx.mesh.create_interval(MPI.COMM_WORLD, nb_cells, [0, l_span])

If you have topologically 1D, but geometrically 2D/3D beams (e.g. beams embedded in 2D/3D space), you will have to use a different method to create mesh. Using gmsh is one way of doing this, but you can also directly handle this with dolfinx. Either way it will not depend on the hermite elements. Something like this should work:

shape = 'interval'
degree = 1
cell = ufl.Cell(shape)
domain= dolfinx.mesh.create_mesh(MPI.COMM_WORLD,
                        cells,
                        x,
                        ufl.Mesh(ufl.VectorElement("Lagrange",cell,degree)) )

There may be a simpler way of doing the above task, but this should be a start.

P.S. - The above code is probably not the best example of a “Minimum Working Example” as there are functions in you module “tools” that other users do not have. For the future, to have the maximum clarity and make debugging easier, please try to strip out the extraneous code and produce a truly minimal example. It may even help your debugging process.

1 Like
3

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!

4

You can use dolfinx’s dofmap ability to extract the relevant dofs and their values. You shouldn’t need to use the rotation to get the x and y displacement, as the rotation is by definition the derivative of the displacement in this formulation. You just have to be careful about the way the dofs are ordered since the rotations are included in the dolfinx/petsc vector. However, they should be pretty easy to tell apart as the rotations will likely be many orders of magnitude lower values that the displacements.