# Creating a 1D beam (full cylinder) clamped-free

Hi all!
I’m struggling to modify this example: Implementation — FEniCSx tutorial
so that I can create a FE model of 1D beam with a cylindrical section (I’m trying to model a wire) clamped on one side and free in the other side.
I guess it’s a fairly simple example when used to fenicsx and dolfinx, but I’m new to it so I would be grateful if someone could give me a basis to work with. Thanks!
Here’s my code right now (not working):

``````from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

L = 1
W = 0.2
mu = 1
rho = 1
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
lambda_ = beta
g = gamma

# nb_cells = 20

gdim, shape, degree = 1, "interval", 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
x = np.array([[0], [0.2], [0.6], [1], [1.3]])
cells = np.array([[0, 1], [1, 2], [2, 3], [3, 4]], dtype=np.int64)
domain = mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain)

V = fem.FunctionSpace(domain, ("Lagrange", 1))

def clamped_boundary(x):
return np.isclose(x[1], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array((0.,) * 1, dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

# u_zero = np.array((0,) * 1, dtype=default_scalar_type)
# print(fem.locate_dofs_geometrical(V, clamped_boundary))
# bc = fem.dirichletbc(u_zero, fem.locate_dofs_geometrical(V, clamped_boundary), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

ds = ufl.Measure("ds", domain=domain)

def epsilon(u):

def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# print(type(uh))
# print(uh.x.array)

with io.XDMFFile(domain.comm, "deformation.xdmf", "w",encoding=io.XDMFFile.Encoding.ASCII) as xdmf:
xdmf.write_mesh(domain)
uh.name = "Deformation"
xdmf.write_function(uh)
``````

This does not make sense.
The input to `clamped boundary` is a `(3, num_points)` array describing points in the xyz plane. You are currently checking if the `y`-value of any point is equal to 0, which is is for all vertices in your mesh (as your mesh is along the x-axis.

Change this to `x[0]`.

There are plenty of other things that doesn’t make sense in your code as well.
Neither `T` nor `f` should be 3D vectors, as you are solving a 1D problem with a scalar variable.

I think you would want to use
`V = fem.FunctionSpace(domain, ("Lagrange", 1, (gdim, )))` to at least get a vector quantity (a 1D vector) for your variational form.

Thank you for your answer. Do you know if some examples of this problem exist? Or something similar. I can’t find anything relevant whereas it seems like a fairly classical example for beginners

You could check: Welcome — Computational Mechanics Numerical Tours with FEniCSx by @bleyerj
which has many examples relating to solid mechanics.

The DOLFINx tutorials (that you refer to above) cover a wide range of PDEs and applications.
There are also several other tutorials:
https://jsdokken.com/fenics22-tutorial/intro.html
and the DOLFINx demos:
https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos.html
covering other aspects.

After some changes I realized that what I wanted was a 1D mesh in a 3D model (forces, displacements, … in 3D). So I get the following code:

``````from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

L = 1
W = 0.2
mu = 1
rho = 1
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
lambda_ = beta
g = gamma

nb_cells = 20

domain = mesh.create_interval(comm=MPI.COMM_WORLD, points=(0, 1), nx=nb_cells)

element_interpolation_order = 1
print(domain.ufl_cell())
V = fem.FunctionSpace(domain, ufl.VectorElement("Lagrange", domain.ufl_cell(), element_interpolation_order, dim=3))

def clamped_boundary(x):
return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array((0.,) * 3, dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

coords = domain.geometry.x

ds = ufl.Measure("ds", domain=domain)

def epsilon(u):

def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
print(u.ufl_shape)
print(V.__dict__)
test = epsilon(u)
print(test)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# print(type(uh))
# print(uh.x.array)

with io.XDMFFile(domain.comm, "deformation.xdmf", "w",encoding=io.XDMFFile.Encoding.ASCII) as xdmf:
xdmf.write_mesh(domain)
uh.name = "Deformation"
xdmf.write_function(uh)
``````

And I have a problem with the shapes of u and v that are (3,) when I run u.ufl_shape. Thus, i have ufl.grad(u) of shape (3,1) and so I get this error in the epsilon function:
ValueError: Cannot take symmetric part of rectangular matrix with dimensions (3, 1)

Any thoughts on this? Thanks!

Well, you really cannot do a symmetric gradient of u if your mesh is not embedded into 3D. Thus, you would have to make your mesh by hand (embed it in 3D):

``````gdim, shape, degree = 3, "interval", 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
x = np.array([[0,0,0], [0.2,0,0], [0.6,0,0], [1,0,0], [1.3,0,0]])
cells = np.array([[0, 1], [1, 2], [2, 3], [3, 4]], dtype=np.int64)
domain = mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain)
``````

and use such a mesh (which now has geometrical dimension 3, check with `print(domain.geometry.dim)`.
Then the gradient is 3D and you get a square matrix for grad(u)

Thanks a lot, the code now seems to work.
Great support, have a nice day!

Just in case it might be useful to someone, here is the code:

``````from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

L = 1
W = 0.2
mu = 1
rho = 1
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
lambda_ = beta
g = gamma

nb_cells = 20

gdim, shape, degree = 3, "interval", 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
x = np.array([[0,0,0], [0.2,0,0], [0.6,0,0], [1,0,0], [1.3,0,0]])
cells = np.array([[0, 1], [1, 2], [2, 3], [3, 4]], dtype=np.int64)
domain = mesh.create_mesh(MPI.COMM_WORLD, cells, x, domain)

# domain = mesh.create_interval(comm=MPI.COMM_WORLD, points=(0, 1), nx=nb_cells)

element_interpolation_order = 1
print(domain.ufl_cell())
V = fem.FunctionSpace(domain, ufl.VectorElement("Lagrange", domain.ufl_cell(), element_interpolation_order, dim=3))

def clamped_boundary(x):
return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array((0.,) * 3, dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

coords = domain.geometry.x

ds = ufl.Measure("ds", domain=domain)

def epsilon(u):

def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
print(V.__dict__)
test = epsilon(u)
print(test)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# print(type(uh))
# print(uh.x.array)

with io.XDMFFile(domain.comm, "deformation.xdmf", "w", encoding=io.XDMFFile.Encoding.ASCII) as xdmf:
xdmf.write_mesh(domain)
uh.name = "Deformation"
xdmf.write_function(uh)
``````

Hey @Paul_Morel (and @dokken)!

First of all, thank you so much for giving the full code since your question is very relevant for my thesis project. I do have a question, when I try to run your code I get the following error: TypeError: Cell._init _() got an unexpected keyword argument ‘geometric dimension’.

I’m very new to programming in general so very sorry if I don’t provide enough information. I tried this on docker with dolfinx:stable and on a virtual machine with Ubuntu 24.04 and I run my (real) machine on Windows 10.

edit: for the VM I simply followed the steps provided on the site to install:

sudo apt update
sudo apt install fenicsx

Here is a version compatible with dolfinx v 0.8.0:

``````from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import basix.ufl

L = 1
W = 0.2
mu = 1
rho = 1
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
lambda_ = beta
g = gamma

nb_cells = 20

gdim, shape, degree = 3, "interval", 1
c_el = basix.ufl.element("Lagrange", shape, degree, shape=(gdim,))
x = np.array([[0, 0, 0], [0.2, 0, 0], [0.6, 0, 0], [1, 0, 0], [1.3, 0, 0]])
cells = np.array([[0, 1], [1, 2], [2, 3], [3, 4]], dtype=np.int64)
domain = mesh.create_mesh(MPI.COMM_WORLD, cells, x, ufl.Mesh(c_el))

# domain = mesh.create_interval(comm=MPI.COMM_WORLD, points=(0, 1), nx=nb_cells)

element_interpolation_order = 1
print(domain.ufl_cell())
V = fem.functionspace(domain, ("Lagrange", element_interpolation_order, (gdim,)))

def clamped_boundary(x):
return np.isclose(x[0], 0)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

u_D = np.array((0.0,) * 3, dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

coords = domain.geometry.x

ds = ufl.Measure("ds", domain=domain)

def epsilon(u):
return ufl.sym(

def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
print(V.__dict__)
test = epsilon(u)
print(test)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

problem = LinearProblem(
a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
uh = problem.solve()

# print(type(uh))
# print(uh.x.array)

with io.XDMFFile(
domain.comm, "deformation.xdmf", "w", encoding=io.XDMFFile.Encoding.ASCII
) as xdmf:
xdmf.write_mesh(domain)
uh.name = "Deformation"
xdmf.write_function(uh)

``````

Thank you so much for taking the time to help me out! This community is great!