Linear elasticity with locally varying constant eigenstrain

variational formulation
1

Hello everyone,

I was searching this board and the internet in general, but it seems I can’t find an answer to my problem. Perhaps because I am also new to FEniCS.
The basic problem is quite simple: I have a periodic 3D box with regularly aligned nodes, and on the basis of linear elasticity I want to calculate the displacements in this box. So far so good, for linear elasticity their are plenty examples and also for periodic boundaries I could find some related topics.

My actual problem is related to the load I want to apply. This is actually a locally varying eigenstrain field, which I read in from a HDF5 file. So I have basically an eigenstrain tensor (3x3) for every DOF, and simply from the theory I would try something like the following:

def eps_el(v, mesh_index):
    eps_tot           = sym(grad(v))
    eigenstrain_local = eigenstrain_field[mesh_index]
    return eps_tot - eigenstrain_local

and then use it e.g. to compute the local stress and so on.
But I do not find a way get this done. I’ve tried to work with Expressions or a dofmap and similar, but perhaps it is just because I am pretty new to the whole FEniCS way of thinking.

Thanks alot for any help!

2

Please make a minimal working code example following the guidelines Read before posting: How do I get my question answered?
Since your data is in a h5 format, i suggest simply doing a 2x2 unit-square mesh, and write the data as a list of numpy nd arrays to illustrate what you would like to do.

Please note that you need to create a mapping from the data you are reading in to the corresponding degree of freedom, as there is no guarantee that your input data is ordered in the same way as the dofs (especially if you want your code to work in parallel).

Also, you will need to use a user-expression to do this efficiently.

3

Thanks for your answer so far, it tells me I do have the right ideas, but can’t figure out how to make them work.

So here is what I did so far:
A) I fill a list with the eigenstrain tensors (here of course only a dummy filling, actually I would map the data read-in from a HDF5 file). Here, at one node I set a tensor unequal to 0.0. The ordering in the list is due to the ordering of the coordinates (loop over x, y, z) in the mesh.

dummy_eigenstrain = np.ndarray(shape=(3, 3), dtype=float)
dummy_eigenstrain.fill(0.0)
dummy_eigenstrain[0][0] = -0.01
eigenstrains = []
for v in vertices(mesh):
    x = v.point().x()
    y = v.point().y()
    z = v.point().z()
    if (x == a0 and y == a0 and z == a0):
        eigenstrains.append(dummy_eigenstrain)
    else:
        eigenstrains.append(0.0*dummy_eigenstrain)

B) I create a UserExpression, which should access the right index in the list:

class getEigenstrain(UserExpression):
    def __init__(self, dx, **kwargs):
        self.dx = dx

    def eval(self, value, 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
        value = eigenstrains[idx]

    def value_shape(self):
        return (3, 3)

eigenstrain_loc = getEigenstrain(dx = a0/2)

C) I try to use this in the calculation of the elastic strain:

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T) - eigenstrain_loc

This gives e.g. the error message:
AttributeError: 'getEigenstrain' object has no attribute '_ufl_shape'

And from here I am pretty much stuck, as everything I try seems to be just random without really understanding what I do.

Anyway, the complete code (still quite long unfortunately) is at the end of this message.

PS.:
Also, as I wrote in the initial message, I am using periodic boundary conditions. I have the feeling, my implementation is not totally right yet (especially at the corners of the domain a test case I did looked quite odd, test case is commented in the code below in the epsilon function). Also, of course, any other hints regarding the code in general are also very welcome.

from __future__ import print_function
from dolfin import *
import numpy as np

# Geometrical and physical parameters
a0      = 2.865
pointsx = 4
pointsy = 4
pointsz = 4
Lx      = a0/2 * pointsx
Ly      = a0/2 * pointsy
Lz      = a0/2 * pointsz

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(pointsx)
Ny      = int(pointsy)
Nz      = int(pointsz)

# Create output file
outputfile = XDMFFile("output.xdmf")

# Create mesh based on hexahedrons
mesh = BoxMesh.create([pmin, pmax], [Nx, Ny, Nz], 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
eigenstrains = []
for v in vertices(mesh):
    x = v.point().x()
    y = v.point().y()
    z = v.point().z()
    if (x == a0 and y == a0 and z == a0):
        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
        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]
        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
        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
        elif near(x[0], xmax): # x+ face
            y[0] = x[0] - xmax
            y[1] = x[1]
            y[2] = x[2]
        elif near(x[1], ymax): # y+ face
            y[0] = x[0]
            y[1] = x[1] - ymax
            y[2] = x[2]
        elif near(x[2], zmax): # z+ face
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - zmax
        else:
            y[0] = -10000
            y[1] = -10000
            y[2] = -10000

# An example what I tried to do with a 'UserExpression'
class getEigenstrain(UserExpression):
    def __init__(self, dx, **kwargs):
        self.dx = dx

    def eval(self, value, 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
        value = eigenstrains[idx]

    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 # - 0.01*Identity(d) #< this was just to test the code
    #return 0.5*(nabla_grad(u) + nabla_grad(u).T) - 0.01*Identity(d) #< this was just to test the code

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(), 1)
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", 1)  
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)
4

I’ve previously shown how to make a tensor-based expression here: User-defined expression for tensor coefficient - #2 by dokken
This should give you some guidelines.

5

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:

fenics_inclusion
fenics_inclusion1266×885 66.1 KB

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)