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

Paul_Morel · February 19, 2024, 3:49pm

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):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


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)

dokken · February 19, 2024, 4:00pm

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.

Paul_Morel · February 19, 2024, 4:08pm

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

dokken · February 19, 2024, 4:12pm

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
http://jsdokken.com/FEniCS23-tutorial/README.html
and the DOLFINx demos:
https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos.html
covering other aspects.

Paul_Morel · February 20, 2024, 10:56am

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):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

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!

dokken · February 20, 2024, 12:23pm

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)

Paul_Morel · February 20, 2024, 2:47pm

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

Paul_Morel · February 20, 2024, 2:49pm

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):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

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)

Dorantor · June 8, 2024, 8:26pm

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 add-apt-repository ppa:fenics-packages/fenics
sudo apt update
sudo apt install fenicsx

dokken · June 9, 2024, 11:39am

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(
        ufl.grad(u)
    )  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)


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)

Dorantor · June 9, 2024, 9:48pm

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

Topic		Replies	Views
Problem to solve a Navier-Stokes problem with a Custom Newton solver in FEniCSx General	6	690	May 23, 2023
Hyperelastic beam dolfinx	5	184	May 28, 2024
Unexpected Behaviour with .sub in fenicsx dolfinx	2	387	September 30, 2022
Magnetostatics example in dolfinx dolfinx	1	2277	May 29, 2020
Solution Error with Lagrange 1 elements dolfinx	2	24	August 16, 2024

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

Related topics