How to add eigenstrain state to linear elasticity

General
1

Hi everyone,

I’m trying to apply an eigenstrain field in a linear elasticity setup and compute the resulting displacements (using the linear thermoelasticity tutorial as a reference). If I understand correctly, I just need to substitute the thermal strain used in the tutorial with a known strain field.

The eigenstrain should be applied only inside a specific region of the mesh. I tested two approaches for defining the eigenstrain: using a DG0 and P1 spaces (options 1 and 2 in the code below). Both run without errors, but the resulting displacement fields differ slightly.

Is one of these approaches considered the recommended way to apply an eigenstrain field? Or is there a more robust/better-practice method for assigning eigenstrains?

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


L, H, W = 5, 0.3, 1
Nx, Ny, Nz = 20, 5, 10
domain = mesh.create_box(
    MPI.COMM_WORLD,
    [(0.0, 0.0, 0.0), (L, H, W)],
    [Nx, Ny, Nz],
    cell_type=mesh.CellType.tetrahedron,
)
gdim = domain.geometry.dim


def lateral_sides(x):
    return np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], L))


E = fem.Constant(domain, 50e3)
nu = fem.Constant(domain, 0.2)
mu = E / 2 / (1 + nu)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)

f = fem.Constant(domain, (0.0, 0.0, 0.0))

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

def sigma(v, eps_star):
    eps_total = eps(v) - eps_star
    return lmbda * ufl.tr(eps_total) * ufl.Identity(gdim) + 2.0 * mu * eps_total
 

# Option 1: define eigenstrain using DG 0
V_strain = fem.functionspace(domain, ("DG", 0, (gdim, gdim)))
eps_star = fem.Function(V_strain, name="eps_star")
num_cells_local = domain.topology.index_map(gdim).size_local
strain_tensor = np.zeros((num_cells_local, gdim, gdim))

def eps_star_mask(x):
    mask = ((1 <= x[0]) & (x[0] <= 4) & (x[1] >= 0.1))
    return mask

strain_cells = dolfinx.mesh.locate_entities(domain, gdim, eps_star_mask)
strain_tensor[strain_cells, 0, 0] = - 0.01 # in x direction
eps_star.x.petsc_vec[:] = strain_tensor.reshape(-1,)


# # Option 2: define eigenstrain using P 1
# V_strain = fem.functionspace(domain, ("P", 1, (gdim, gdim)))
# eps_star = fem.Function(V_strain, name="eps_star")
# dof_coords = V_strain.tabulate_dof_coordinates()
# x_coords = dof_coords[:, 0]
# y_coords = dof_coords[:, 1]
# z_coords = dof_coords[:, 2]
# mask = ((1 <= x_coords) & (x_coords <= 4) & (y_coords >= 0.1))
# strain_values = np.zeros((len(x_coords), gdim, gdim))
# strain_values[mask, 0, 0] = - 0.01 # in x direction
# eps_star.x.petsc_vec[:] = strain_values.reshape(-1,)


Vu = fem.functionspace(domain, ("Lagrange", 1, (gdim,)))
du = ufl.TrialFunction(Vu)
u_ = ufl.TestFunction(Vu)
Wint = ufl.inner(sigma(du, eps_star), eps(u_)) * ufl.dx
aM = ufl.lhs(Wint)
LM = ufl.rhs(Wint) + ufl.inner(f, u_) * ufl.dx

lateral_dofs_u = fem.locate_dofs_geometrical(Vu, lateral_sides)
bcu = [fem.dirichletbc(np.zeros((gdim,)), lateral_dofs_u, Vu)]

u = fem.Function(Vu, name="Displacement")
problem = fem.petsc.LinearProblem(aM, LM, u=u, bcs=bcu)
problem.solve()

print(f"max displacement: {np.round(np.max(np.abs(u.x.array)), 5)}")

with io.XDMFFile(domain.comm, "strain.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(eps_star)

with io.XDMFFile(domain.comm, "displacement.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(u)
2

Technically, it looks to me like the eigenstrain lives in a tensor L^2 space, hinting at DG0 being the most appropriate space. This would also allow you to localize it to subdomains. But H^1 \subset L^2 so I don’t think there is a fundamental issue with choosing CG1.

Importantly, you’re not solving for this field. The choice of function space typically becomes precarious for the test and trial spaces. For the PDE parameters such as material coefficients (which this eigenstrain effectively is) the choice is far less crucial. Typically, it only affects “eyeball norm” quality on coarse meshes.

1 Like
3

Hi @Stein

thanks a lot for the clear explanation.

I was confused because in the expression eps_total = eps(v) - eps_star, the eps_star function has different array sizes (eps_star.x.array.shape) depending on the function space:

P1 → shape (12474,)
DG0 → shape (54000,)

But I assume this is not a problem and FEniCSx knows how to correctly interpret both in the assembly process?

4

Yes absolutely. As a form fenics thinks of the objects as functions (which can be added, multiplied, …) and not as their arrays of basisfunction coefficients.

1 Like
5

Just go with the regularity/continuity of the data as you have it available. If this is from a temperature change, as in Jeremy’s demo, then the thermal strains are in continuous Lagrange space. But if you have the thermal expansion coefficient discontinuous, because there are subdomains of different materials in the domain, then make the whole thermal strain tensor discontinuous. Different strains have different regularities based on their physical origin.

FEniCSx will evaluate all objects in the dolfinx.fem.form at quadrature points, so you do not need to care about them having different discretization. Note that the lowest regularity you can input here would be a Quadrature element (not the DG0 space), which is often the right choice for e.g. plastic strains.

2 Likes
6

Hi Michal,

thanks a lot for the clarification :+1: