Thank you, @dokken! This really solves a big part of my problem, but something is still not working.

I have changed my expression and strain to

```
# Acces the eigenstrains via UserExpression
class getEigenstrain(UserExpression):
def __init__(self, dx, **kwargs):
# Call superclass constructor with keyword arguments
super().__init__(**kwargs)
self.dx = dx
def eval(self, values, x):
px = int( x[0] / self.dx )
py = int( x[1] / self.dx )
pz = int( x[2] / self.dx )
idx = px + py * Nx + pz * Nx*Ny
eigenstrain_loc = eigenstrains[idx]
values[0] = eigenstrain_loc[0][0]
values[1] = eigenstrain_loc[0][1]
values[2] = eigenstrain_loc[0][2]
values[3] = eigenstrain_loc[1][0]
values[4] = eigenstrain_loc[1][1]
values[5] = eigenstrain_loc[1][2]
values[6] = eigenstrain_loc[2][0]
values[7] = eigenstrain_loc[2][1]
values[8] = eigenstrain_loc[2][2]
def value_shape(self):
return (3, 3)
eigenstrain_loc = getEigenstrain(dx = a0/2)
# Define strain and stress
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T) - eigenstrain_loc
```

which works great.

But after a few tests I realised, that my results do not look like I expected them.

For example:

I place a 3x3x3 eigenstrain “inclusion” centred into a 5x5x5 domain, and the eigenstrain tensor mimics some kind of volumetric expansion (eps11 = eps22 = ep33 > 0). And what I get is this:

As you may see, the displacement vectors coming from the “inclusion” only point into negative coordinate directions. I tested this in different variations (domain size, inclusion size, elements, …) and it is always the same picture.

So I was thinking that it either has to do with my periodic boundary implementation + Lagrange multiplier (for which I mainly followed https://comet-fenics.readthedocs.io/en/latest/demo/periodic_homog_elas/periodic_homog_elas.html) or the expression does not fully work as expected.

My updated code:

```
from __future__ import print_function
from dolfin import *
import numpy as np
# Geometrical and physical parameters
a0 = 2.865
cellsx = 4
cellsy = 4
cellsz = 4
Lx = a0/2 * cellsx
Ly = a0/2 * cellsy
Lz = a0/2 * cellsz
E = 210e3
nu = 0.3
xmin = 0.0
ymin = 0.0
zmin = 0.0
pmin = Point(float(xmin), float(ymin), float(zmin))
xmax = Lx
ymax = Ly
zmax = Lz
pmax = Point(float(xmax), float(ymax), float(zmax))
Nx = int(cellsx) + 1
Ny = int(cellsy) + 1
Nz = int(cellsz) + 1
# Create output file
outputfile = XDMFFile("output.xdmf")
# Create mesh based on hexahedrons
mesh = BoxMesh.create([pmin, pmax], [Nx-1, Ny-1, Nz-1], CellType.Type.hexahedron)
outputfile.write(mesh)
# Fill eigenstrains due to mesh indices
# Actually here I would have read-in a HDF5 file, ordered by xyz,
# and would try to map the values to the respective DOFs
dummy_eigenstrain = np.ndarray(shape=(3, 3), dtype=float)
dummy_eigenstrain.fill(0.0)
dummy_eigenstrain[0][0] = +0.01
dummy_eigenstrain[1][1] = +0.01
dummy_eigenstrain[2][2] = +0.01
eigenstrains = []
for v in vertices(mesh):
x = v.point().x()
y = v.point().y()
z = v.point().z()
if (x >= a0/2 and x <= 3*a0/2 and y >= a0/2 and y <= 3*a0/2 and z >= a0/2 and z <= 3*a0/2):
eigenstrains.append(dummy_eigenstrain)
else:
eigenstrains.append(0.0*dummy_eigenstrain)
# Define periodic boundary condition
class PeriodicBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool( ( near(x[0], xmin) or near(x[1], ymin) or near(x[2], zmin) )
and
( not ((near(x[2], zmin) and near(x[0], xmax)) or
(near(x[1], ymin) and near(x[0], xmax)) or
(near(x[0], xmin) and near(x[2], zmax)) or
(near(x[1], ymin) and near(x[2], zmax)) or
(near(x[0], xmin) and near(x[1], ymax)) or
(near(x[2], zmin) and near(x[1], ymax)) ))
and on_boundary )
def map(self, x, y):
if near(x[0], xmax) and near(x[1], ymax) and near(x[2], zmax): # x+, y+, z+ point
y[0] = x[0] - xmax
y[1] = x[1] - ymax
y[2] = x[2] - zmax
#print("%i %i %i --> %i %i %i" % (x[0] / (a0/2), x[1] / (a0/2), x[2] / (a0/2), y[0] / (a0/2), y[1] / (a0/2), y[2] / (a0/2)))
elif near(x[0], xmax) and near(x[1], ymax): # x+, y+ edge
y[0] = x[0] - xmax
y[1] = x[1] - ymax
y[2] = x[2]
#print("%i %i %i --> %i %i %i" % (x[0] / (a0/2), x[1] / (a0/2), x[2] / (a0/2), y[0] / (a0/2), y[1] / (a0/2), y[2] / (a0/2)))
elif near(x[0], xmax) and near(x[2], zmax): # x+, z+ edge
y[0] = x[0] - xmax
y[1] = x[1]
y[2] = x[2] - zmax
#print("%i %i %i --> %i %i %i" % (x[0] / (a0/2), x[1] / (a0/2), x[2] / (a0/2), y[0] / (a0/2), y[1] / (a0/2), y[2] / (a0/2)))
elif near(x[1], ymax) and near(x[2], zmax): # y+, z+ edge
y[0] = x[0]
y[1] = x[1] - ymax
y[2] = x[2] - zmax
#print("%i %i %i --> %i %i %i" % (x[0] / (a0/2), x[1] / (a0/2), x[2] / (a0/2), y[0] / (a0/2), y[1] / (a0/2), y[2] / (a0/2)))
elif near(x[0], xmax): # x+ face
y[0] = x[0] - xmax
y[1] = x[1]
y[2] = x[2]
#print("%i %i %i --> %i %i %i" % (x[0] / (a0/2), x[1] / (a0/2), x[2] / (a0/2), y[0] / (a0/2), y[1] / (a0/2), y[2] / (a0/2)))
elif near(x[1], ymax): # y+ face
y[0] = x[0]
y[1] = x[1] - ymax
y[2] = x[2]
#print("%i %i %i --> %i %i %i" % (x[0] / (a0/2), x[1] / (a0/2), x[2] / (a0/2), y[0] / (a0/2), y[1] / (a0/2), y[2] / (a0/2)))
elif near(x[2], zmax): # z+ face
y[0] = x[0]
y[1] = x[1]
y[2] = x[2] - zmax
#print("%i %i %i --> %i %i %i" % (x[0] / (a0/2), x[1] / (a0/2), x[2] / (a0/2), y[0] / (a0/2), y[1] / (a0/2), y[2] / (a0/2)))
else:
y[0] = x[0]
y[1] = x[1]
y[2] = x[2]
# Acces the eigenstrains via UserExpression
class getEigenstrain(UserExpression):
def __init__(self, dx, **kwargs):
# Call superclass constructor with keyword arguments
super().__init__(**kwargs)
self.dx = dx
def eval(self, values, x):
px = int( x[0] / self.dx )
py = int( x[1] / self.dx )
pz = int( x[2] / self.dx )
idx = px + py * Nx + pz * Nx*Ny
eigenstrain_loc = eigenstrains[idx]
values[0] = eigenstrain_loc[0][0]
values[1] = eigenstrain_loc[0][1]
values[2] = eigenstrain_loc[0][2]
values[3] = eigenstrain_loc[1][0]
values[4] = eigenstrain_loc[1][1]
values[5] = eigenstrain_loc[1][2]
values[6] = eigenstrain_loc[2][0]
values[7] = eigenstrain_loc[2][1]
values[8] = eigenstrain_loc[2][2]
#if values[0] > 0.0:
#print("%i %i %i LOCAL tensor\n" % (px, py, pz), values.reshape(self.value_shape()))
def value_shape(self):
return (3, 3)
eigenstrain_loc = getEigenstrain(dx = a0/2)
# Define strain and stress
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T) - eigenstrain_loc
def sigma(u):
lmbda = E*nu/((1.0+nu) * (1.0-2.0*nu))
mu = E/(2.0 * (1.0+nu))
return lmbda*tr(epsilon(u))*Identity(d) + 2.0*mu*(epsilon(u))
VE = VectorElement("CG", mesh.ufl_cell(), 2)
RE = VectorElement("R", mesh.ufl_cell(), 0)
VLFS = FunctionSpace(mesh, MixedElement([VE, RE]), constrained_domain=PeriodicBoundary())
VFS = FunctionSpace(mesh, VE)
# Define variational problem
u, lamb = TrialFunctions(VLFS)
d = u.geometric_dimension() # space dimension
v, lamb_ = TestFunctions(VLFS)
x_res = Function(VLFS)
F = inner(sigma(u), epsilon(v)) * dx
a, L = lhs(F), rhs(F)
a += dot(lamb_, u)*dx + dot(lamb, v)*dx
# Compute solution
#solve(a == L, x_res, [])
solve(a == L, x_res, [], solver_parameters={"linear_solver": "cg"})
u, lamb = split(x_res)
# Create function space
FS = FunctionSpace(mesh, "CG", 2)
u = project(u, VFS)
u.rename("u", "displacement vector")
# Write solution to xdmf
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, FS)
outputfile.write(u, 0.0)
outputfile.write(u_magnitude, 0.0)
```