How to parameterize the linear elasticity example from tutorial?

5 months ago by
Hello all,

I am new to FeNiCS and beginner at coding in general. I was able to follow the Linear Elasticity example from the tutorial (3.3), and now I am trying to modify the code such that I can model the displacement according to variable values of young's modulus/poisson ratio. I am new to coding, so I am not sure how I can begin to change my constants into parameterized variables. Could someone show me an example program where this is done, or give me a push in the right direction?

Thank you,
Community: FEniCS Project
Do you want the material parameters to vary in space or dependent on the strain/displacement or do you want something completely different?
written 5 months ago by klunkean  
The second one, I would like the material parameters to be dependent on the displacement in the z direction.
written 5 months ago by Joshua Lee  

1 Answer

5 months ago by
If you can express your material parameter as a closed mathematical expression, you can just use python functions. Note, however, that the problem will become non-linear. Consider the following MWE:
from fenics import *

mesh = UnitCubeMesh(4,4,4)
FE = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
VFE = VectorElement(FE)
V = FunctionSpace(mesh, VFE)

u = Function(V)
delu = TrialFunction(V)
v = TestFunction(V)

def bottom(x, on_boundary):
    return near(x[2],0.)
def top(x, on_boundary):
    return near(x[2],1.)

dbc = [DirichletBC(V, Constant((0.,0.,0.)), bottom), 
       DirichletBC(V.sub(2), Constant(.1), top)]

eye = Identity(3)
eyOey = outer(eye,eye)
eye4 = as_tensor(.5*(eye[i,k]*eye[j,l]+eye[i,l]*eye[j,k]), (i,j,k,l))
D4 = eye4 - 1./3.*eyOey

# bulk and shear modulus, shear modulus as function of u
K = 160e+9

def G(u):   
    # just an example
    G0 = 80e+9
    uzMax = .1
    return G0*(1.-u[2]/uzMax)

#stiffness tensor
def C4(u):
    return K*eyOey + G(u)*D4

eps = sym(grad(u))
def sigma(u):
    return as_tensor(C4(u)[i,j,k,l]*eps[k,l], (i,j))

#nonlinear form
F = inner(grad(v),sigma(u))*dx

#gateaux derivative
J = derivative(F, u, delu)

solve(F==0, u, dbc, J=J)

File('disp.pvd') << u​
Thank you! This is very helpful. Could you please explain what the

"eye = Identity(3)
eyOey = outer(eye,eye)
eye4 = as_tensor(.5*(eye[i,k]*eye[j,l]+eye[i,l]*eye[j,k]), (i,j,k,l))
D4 = eye4 - 1./3.*eyOey"

part of the code is doing?

I also tried running this code through Docker, but I got the error message:

"SyntaxError: Non-ASCII character '\xe2' in file on line 52, but no encoding declared; see for details"

I tried following the link, but I just got sent to a 404. Do you know why I got this error?

Thank you!

written 5 months ago by Joshua Lee  
The four lines define the following quantities:
eye - Second order unit tensor  $\mathbf{1}$1
eyOey - Dyadic product of second order unit tensor with itself $\mathbf{1}\otimes\mathbf{1}$11  (Trace projector)
eye4 -  Fourth order symmetric unit tensor $\mathbf{1}^{4,\mathrm{s}}=\frac{1}{2}\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)\mathbf{e}_i\otimes\mathbf{e}_j\otimes\mathbf{e}_k\otimes\mathbf{e}_l$14,s=12 (δikδjl+δilδjk)eiejekel
D4 - Fourth order deviatoric projector  $\mathbf{D}^{4,\mathrm{s}}=\mathbf{1}^{4,\mathrm{s}}-\frac{1}{3}\mathbf{1}\otimes\mathbf{1}$D4,s=14,s13 11

With these quantities we can write Hooke's tangent tensor as

Where  $K$K  and  $G$G  are bulk and shear modulus respectively.

As for your error: This is a simple copy-paste error. I don't know what happens in detail, but try to go to the very last line of the pasted code and delete the very last character. There is some invisible character which gets pasted and does not belong there.
written 5 months ago by klunkean  
Hi klunkean,

Could you help me with another issue? I am having trouble extracting the sigma function, similar as to how the original tutorial plotted the stress intensity. I would like to plot the stress intensity as a function of displacement.

Thank you for your help!
written 4 months ago by Joshua Lee  
Hi again.

If your stress is declared as above you could project it onto a finite element space after computing the displacements. Like this:
solve(F==0, u, dbc, J=J)
tensorFE = TensorElement(FiniteElement('DG', mesh.ufl_cell(), 0))
tensorFunSpace = FunctionSpace(mesh,tensorFE)

stress = project(sigma(u), tensorFunSpace)
stressFile = File('stress.pvd')
stressFile << stress​

Then you can view stress.pvd in paraview. Hope that helps.
written 4 months ago by klunkean  
you have been very helpful, thank you again for your help!
written 4 months ago by Joshua Lee  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »