2D linear elasticity using dolfinx

General
1

Hello everyone,
I’m currently redo the example 2D linear elasticity from this link “2D linear elasticity — Numerical tours of continuum mechanics using FEniCS master documentation” in which they use dolfin.

When i try to rewrite the code in dolfinx, i got the error as below:
Symmetric part of tensor with rank !=2 is undefined.
My code:
from mpi4py import MPI
from dolfinx import mesh, fem, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import sym, grad, tr, Identity, TrialFunction, TestFunction, inner, dx

L, H, Nx, Ny = 25., 1., 250, 10
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[[0.0, 0.0], [L, H]],
[Nx, Ny],
mesh.CellType.quadrilateral,
)

def eps(v):
return sym(grad(v))

E = fem.Constant(domain, 1e5)
nu = fem.Constant(domain, 0.3)
model = “plane_stress”

mu = E/2/(1+nu)
lmbda = Enu/(1+nu)/(1-2nu)
if model ==“plane_stress”:
lmbda = 2mulmbda/(lmbda+2*mu)

def sigma(v):
return lmbda*tr(eps(v))Identity(2) + 2.0 * muesp(v)

rho_g = 1e-3
f = fem.Constant(domain, (0, -rho_g))
V = fem.functionspace(domain, (“Lagrange”, 2))

du= TrialFunction(V)
u_= TestFunction(V)
a = inner (sigma(du), eps(u_)) *dx
L = inner(f, u_)*dx

def left(x, on_boundary):
return x[0] < fem.DOLFINX_EPS

bc = fem.DirichletBC(V, fem.Constant((0.,0.)), left)
u = fem.Function(V, name=“Displacement”)

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

plot(1e3uh, mode=“displacement”)
print(“Maximal deflection:”, -uh(L,H/2.)[1])
print(“Beam theory deflection:”, float(3rho_g*L4/2/E/H3))

If you experienced or had the solution, please share me.
Thank you so much for your reading!

2

Please post the code properly, i.e. within a block delimited by “```”.

Your mistake is

V = fem.functionspace(domain, (“Lagrange”, 2))

It should rather be

V = fem.functionspace(domain, (“Lagrange”, 2, (2,)))

where the former defines a FE space for an unknown scalar field, while the latter defines it for a unknown vector field in 2D.
P

1 Like
3

I’m so sorry about my mistake.

Thank you so much for your help (it worked).

4

Good morning
Trying to reproduce this code to understand how to retrieve displacements with dolfinx, I have an issue at this line :

bc = fem.DirichletBC(V, fem.Constant((0.,0.)), left)

says :

TypeError: Constant.__init__() missing 1 required positional argument: 'c'

I tried a few things without success :

# bc = fem.DirichletBC(V, fem.Constant([0.0, 0.0]), left)
# bc = fem.dirichletbc(V, fem.Constant((0.,0.)), left)
# bc = dolfinx.DirichletBC(V, fem.Constant((0.,0.)), left)
# constant_zero = fem.Constant(domain, 0.0)
# bc = fem.DirichletBC(V, constant_zero, left)

Can you please guide me to the right syntax ?
Here bellow is the full code :

import os
arch = os.getenv("ARGS", "real")

try:
    import google.colab  # noqa: F401
except ImportError:
    import ufl
    import dolfinx
else:
    try:
        import ufl
        import dolfinx
    except ImportError:
        if arch != "complex":
            !wget "https://fem-on-colab.github.io/releases/fenicsx-install-real.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
        else:
            !wget "https://fem-on-colab.github.io/releases/fenicsx-install-complex.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
        import ufl
        import dolfinx

from mpi4py import MPI
from dolfinx import mesh, fem, plot
from dolfinx.fem import dirichletbc
from dolfinx.fem.petsc import LinearProblem
from ufl import sym, grad, tr, Identity, TrialFunction, TestFunction, inner, dx

L, H, Nx, Ny = 25., 1., 250, 10
domain = mesh.create_rectangle(MPI.COMM_WORLD,
[[0.0, 0.0], [L, H]],
[Nx, Ny],
mesh.CellType.quadrilateral,
)

def eps(v):
   return sym(grad(v))

E = fem.Constant(domain, 1e5)
nu = fem.Constant(domain, 0.3)
model = "plane_stress"

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model =="plane_stress":
   lmbda = 2*mu*lmbda/(lmbda+2*mu)

def sigma(v):
   return lmbda*tr(eps(v))*Identity(2) + 2.0 * mu*eps(v)

rho_g = 1e-3
f = fem.Constant(domain, (0, -rho_g))
V = fem.functionspace(domain, ("Lagrange", 2, (2,)))

du= TrialFunction(V)
u_= TestFunction(V)
a = inner (sigma(du), eps(u_)) *dx
L = inner(f, u_)*dx

def left(x, on_boundary):
   return x[0] < fem.DOLFINX_EPS

bc = fem.DirichletBC(V, fem.Constant((0.,0.)), left)
# bc = fem.DirichletBC(V, fem.Constant([0.0, 0.0]), left)
# bc = fem.dirichletbc(V, fem.Constant((0.,0.)), left)
# bc = dolfinx.DirichletBC(V, fem.Constant((0.,0.)), left)
# constant_zero = fem.Constant(domain, 0.0)
# bc = fem.DirichletBC(V, constant_zero, left)

u = fem.Function(V, name="Displacement")

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

plot(1e3*uh, mode="displacement")
print("Maximal deflection:", -uh(L,H/2.)[1])
print("Beam theory deflection:", float(3*rho_g*L4/2/E/H3))
5

Look at all the other times the code uses fem.Constant: you’ll see that it expects the mesh as the first argument.

1 Like
6

Good morning and thanks for the tip.
The following seems to work :

from petsc4py import PETSc  # Import PETSc

# Locate left boundary DOFs
left_dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0))

# Create a constant value for the boundary condition
constant_zero = fem.Constant(domain, np.array([0.0, 0.0], dtype=PETSc.ScalarType))

# Create the Dirichlet boundary condition using dirichletbc function
bc = fem.dirichletbc(constant_zero, left_dofs, V)

Is it a correct way to do ?

Then I run into other issues I am currently trying to handle on my own for now…

7

Yes, you correctly defined the constant in dolfinx by including the mesh as the first argument. What issue are you encountering now? Can I have your contact information? I’m working on solid mechanics, which I believe is also your area of interest.

8

Thanks for confirmation

What issue are you encountering now

Well first :

----> plot(1e3*uh, mode="displacement")
TypeError: 'module' object is not callable

…then :

----> 1 print("Maximal deflection:", -uh(L,H/2.)[1])

2 frames

/usr/local/lib/python3.10/dist-packages/ufl/core/terminal.py in evaluate(self, x, mapping, component, index_values, derivatives)
     41     def evaluate(self, x, mapping, component, index_values, derivatives=()):
     42         """Get *self* from *mapping* and return the component asked for."""
---> 43         f = mapping.get(self)
     44         # No mapping, trying to evaluate self as a constant
     45         if f is None:

AttributeError: 'float' object has no attribute 'get'
9

These are not dolfinx functions.
Please consider the DOLFINx tutorial and demos to understand how to plot (using pyvista) and how to evaluate at points
https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#plotting-the-solution-over-a-line
https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain

10

Ok now it works all fine :

Maximal deflection: -3.12e-04
Beam theory : -3.11e-04

Full code goes like this :

import os
arch = os.getenv("ARGS", "real")

try:
    import google.colab  # noqa: F401
except ImportError:
    import ufl
    import dolfinx
else:
    try:
        import ufl
        import dolfinx
    except ImportError:
        if arch != "complex":
            !wget "https://fem-on-colab.github.io/releases/fenicsx-install-real.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
        else:
            !wget "https://fem-on-colab.github.io/releases/fenicsx-install-complex.sh" -O "/tmp/fenicsx-install.sh" && bash "/tmp/fenicsx-install.sh"
        import ufl
        import dolfinx

from mpi4py import MPI
from dolfinx import mesh, fem, plot
from dolfinx.fem import dirichletbc, locate_dofs_geometrical, Function
from dolfinx.fem.petsc import LinearProblem
from ufl import sym, grad, tr, Identity, TrialFunction, TestFunction, inner, dx
import numpy as np                                       # Numerical operations
from petsc4py import PETSc  # Import PETSc
import pyvista

# Define the length, height, and number of divisions of the beam
length, height, Nx, Ny = 12., 1., 24, 2

# Create a mesh of the beam using the number of divisions and the length and height
domain = mesh.create_rectangle(MPI.COMM_WORLD,
                               [[0.0, 0.0], [length, height]],
                               [Nx, Ny],
                               mesh.CellType.quadrilateral,)

# Define the strain-displacement relationship
def eps(v):
   return sym(grad(v))

# Define the Young's modulus and Poisson's ratio of the beam
E_val = 1e5
E = fem.Constant(domain, E_val)
nu = fem.Constant(domain, 0.3)

# Define the type of beam model
model = "plane_stress"

# Calculate the Lame parameters
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

# Adjust the Lame parameters for the plane stress model
if model =="plane_stress":
   lmbda = 2*mu*lmbda/(lmbda+2*mu)

# Define the stress-strain relationship
def sigma(v):
   return lmbda*tr(eps(v))*Identity(2) + 2.0 * mu*eps(v)

# Define the density and the gravitational acceleration
rho_g = 1e-3
f = fem.Constant(domain, (0, -rho_g))

# Define the function space for the displacement
V = fem.functionspace(domain,
                      ("Lagrange",
                       2,                            # polynomial degree
                       (2,)))                        # function space dimension

# Define the trial and test functions for the displacement
du= TrialFunction(V)
u_= TestFunction(V)

# Define the variational formulation for the equilibrium equation
a = inner (sigma(du), eps(u_)) *dx
L = inner(f, u_)*dx

# Locate left boundary DOFs
left_dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.0))

# Create a constant value for the boundary condition
constant_zero = fem.Constant(domain, np.array([0.0, 0.0], dtype=PETSc.ScalarType))

# Create the Dirichlet boundary condition using dirichletbc function
bc = fem.dirichletbc(constant_zero, left_dofs, V)

# Define the function for the displacement
u = fem.Function(V, name="Displacement")

# Define the linear problem for the equilibrium equation
problem = LinearProblem(a,
                        L,
                        bcs=[bc],
                        petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

# Solve the linear problem
uh = problem.solve()

# plot domain
pyvista.start_xvfb()
tdim = domain.topology.dim             # = 2 topological dimensions of the mesh
domain.topology.create_connectivity(tdim, tdim)
topology, cell_types, geometry = plot.vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

plotter = pyvista.Plotter(notebook=True)
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
plotter.show(jupyter_backend="html", return_viewer=True)

# plot uy

u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
uy = uh.x.array.real[1::2]
u_grid.point_data["uy"] = uy
u_grid.set_active_scalars("uy")
u_plotter = pyvista.Plotter(notebook=True)
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
u_plotter.show_axes()
u_plotter.show(jupyter_backend="html", return_viewer=True)

print("Maximal deflection: {:.2e}".format(min(uy)))
beam_theory_deflection = (-rho_g * length**4) / (8*E_val*height**3/12)
print("Beam theory : {:.2e}".format(beam_theory_deflection))