Non-uniform Young's modulus for linear elasticity

General
1

Hi, I tried to replicate the basic linear elasticity code by Jeremy Bleyer in this link with addition of trying to use a non-uniform E as my trial example.

Tried this function as a simple non-uniform E, however still got error as value error of invalid ranks in product when computing sigma. I suspect the single line when defining my E is still incorrect. Any functions that I can use to correct this bug? Thank you for the help.

E = 210E3 + 21E3 * sin (\pi * x)

Now I’m using the E = fem.Constant(domain, 210e3 + 21e3 * np.sin(np.pi * xc)) as the function to define E to result the array with xc also an array of x-coordinates (domain.geometry.x[:,0]).

import numpy as np
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction, nabla_grad

from mpi4py import MPI

from dolfinx import fem, io
import dolfinx.fem.petsc
from dolfinx.mesh import create_rectangle, CellType

length, height = 10, 1.0
Nx, Ny = 50, 5
domain = create_rectangle(
    MPI.COMM_WORLD,
    [np.array([0, 0]), np.array([length, height])],
    [Nx, Ny],
    cell_type=CellType.quadrilateral,
)

dim = domain.topology.dim
print(f"Mesh topology dimension d={dim}.")

degree = 2
shape = (dim,)  # this means we want a vector field of size `dim`
V = fem.functionspace(domain, ("P", degree, shape))

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

xc = domain.geometry.x[:,0]
E = fem.Constant(domain, 210e3 + 21e3 * np.sin(np.pi * xc))
nu = fem.Constant(domain, 0.3)

lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)


def epsilon(v):
    return 0.5*(nabla_grad(v) + nabla_grad(v).T)


def sigma(v):
    return lmbda * tr(epsilon(v)) * Identity(dim) + 2 * mu * epsilon(v)

print("mu (UFL):\n", mu)
print("epsilon (UFL):\n", epsilon(u_sol))
print("sigma (UFL):\n", sigma(u_sol))

I used Fenicsx version 0.8.0 and Python 3.12.2 in a MBP M1.

2 
xc = domain.geometry.x[:,0]
E = fem.Constant(domain, 210e3 + 21e3 * np.sin(np.pi * xc))

This isn’t a valid constant. Consider a UFL expression or interpolating your function into an appropriate space.

3

I tried to use the fem.Function.expressions as a change but found out that the Expression and Trace has an unsupported operand type. Any suggestion on this? Thank you.

x = ufl.SpatialCoordinate(domain)
E = 210e3 + 21e3 * ufl.classes.Sin(x[0] * np.pi)
lmbd = E * nu / (1 + nu) / (1 - 2 * nu)
mu_ = E / 2 / (1 + nu)

lmbda = fem.Expression(lmbd, V.element.interpolation_points())
mu = fem.Expression(mu_, V.element.interpolation_points())

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

def sigma(v):
    return lmbda * tr(epsilon(v)) * Identity(dim) + 2 * mu * epsilon(v)

TypeError: unsupported operand type(s) for *: ‘Expression’ and ‘Trace’

4

If you want to go down the fem.Expression route, you should interpolate lmbda and mu into a FE space before using them in your python function sigma.

5

Now the code is working. Adding a FE space to act as interpolation.

E = 210e3 + 21e3 * ufl.classes.Sin(x[0] * np.pi)
lmbd = E * nu / (1 + nu) / (1 - 2 * nu)
mu_ = E / 2 / (1 + nu)

W = fem.functionspace(domain, ("Discontinuous Lagrange", 0))

lmbda_exp = fem.Expression(lmbd, W.element.interpolation_points())
mu_exp = fem.Expression(mu_, W.element.interpolation_points())

lmbda = fem.Function(W)
mu = fem.Function(W)

lmbda.interpolate(lmbda_exp)
mu.interpolate(mu_exp)

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

def sigma(v):
    return lmbda * tr(epsilon(v)) * Identity(dim) + 2 * mu * epsilon(v)