How to view/ print the stress/ strain tensor

Question

Hi, I want to view my stress tensor and I used the following. However, it keeps giving me a runtime error.

E = 200*(10**9)
nu = 0.24

mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu))

# Stress
sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(mesh.geometry().dim())   

T = TensorFunctionSpace(mesh, 'CG', 1)
stress = project(sigma(u),T)
File("stress.pvd") << stress
M = stress.vector().array()
print M

Should I use 'Discontinuous Lagrange'? Your help is highly appreciated. Thanks.

I checked this link too, but couldn't figure it out.
http://fenicsproject.org/qa/4051/how-to-get-numpy-array-of-a-tensor
http://fenicsproject.org/qa/5569/evaluating-a-tensor-valued-function-at-a-point?show=5569#q5569

Answer 1

Hi,

write
stress = project(sigma,T)

instead of
stress = project(sigma(u),T)

Comments

I tried. It still gives me this error:

Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "trial.py", line 32, in <module>
    stress = project(sigma,T)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/projection.py", line 111, in project
    form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py", line 254, in assemble_system
    assembler = cpp.SystemAssembler(A_dolfin_form, b_dolfin_form, bcs)
  File "/usr/lib/python2.7/dist-packages/dolfin/cpp/fem.py", line 2338, in __init__
    _fem.SystemAssembler_swiginit(self,_fem.new_SystemAssembler(*args))
RuntimeError: 

---

Could you upload the whole code please?

---

This is the code:

from dolfin import *
d0 = 8.
d1 = 1.
d2 = 1.
n0 = int(d0*6)
n1 = int(d1*6)
n2 = int(d2*6)
o = Point(0, 0, 0)
l = Point(8, 1, 1)
Mesh = BoxMesh(o, l, n0, n1, n2)
Mesh.init()

V  = VectorFunctionSpace(Mesh, "CG", 1)

tol = 1E-14 # tolerance 

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 8.0)      

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 1.0)

left = Left()
top = Top()
right = Right()

boundaries = FacetFunction("uint", Mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)      

u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, -1.9, 0.0))

E = 200*(10**9)
nu = 0.24

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

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

def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d) 

T = TensorFunctionSpace(Mesh, 'CG', 1)
stress = project(sigma,T)
File("stress.pvd") << stress
M = stress.vector().array()
print M

ds = Measure("ds")[boundaries]

g_T= Constant((0.0, -8500000.0, 0.0))

a = inner(grad(v),sigma(u))*dx 

L= dot(v,g_T)*ds(2)+dot(v,g_T)*dx

left_end_displacement = Constant((0.0, 0.0, 0.0))
bc = [DirichletBC(V, left_end_displacement, left)] 
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bc)
solver = LinearVariationalSolver(problem)
solver.parameters["symmetric"] = True
solver.solve()
plot(u, mode="displacement",interactive=True)

---

at the end of your code, redefine v, because in functioni propose to change the lines
a = inner(grad(v),sigma(u))*dx
L= dot(v,g_T)ds(2)+dot(v,g_T)dx

with

w = TestFunction(V)
a = inner(grad(w),sigma(u))*dx
L= dot(w,g_T)ds(2)+dot(w,g_T)dx

here i send the whole code:

from dolfin import *
d0 = 8.
d1 = 1.
d2 = 1.
n0 = int(d06)
n1 = int(d16)
n2 = int(d2*6)
o = Point(0, 0, 0)
l = Point(8, 1, 1)
Mesh = BoxMesh(o, l, n0, n1, n2)
Mesh.init()

V = VectorFunctionSpace(Mesh, "CG", 1)

tol = 1E-14 # tolerance

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 8.0)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 1.0)

left = Left()
top = Top()
right = Right()

boundaries = FacetFunction("uint", Mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)

u = TrialFunction(V)
u = Function(V)
f = Constant((0.0, -1.9, 0.0))

E = 200*(10**9)
nu = 0.24

mu, lmbda = Constant(E/(2(1.0 + nu))), Constant(Enu/((1 + nu)(1 - 2nu)))

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

def sigma(v):
def sigma(v):
return 2.0musym(grad(v)) + lmbdatr(sym(grad(v)))Identity(v.cell().topological_dimension() )

T = TensorFunctionSpace(Mesh, 'CG', 1)
stress = project(sigma,T)
File("stress.pvd") << stress
M = stress.vector().array()
print M

ds = Measure("ds")[boundaries]

g_T= Constant((0.0, -8500000.0, 0.0))

w = TestFunction(V)
a = inner(grad(w),sigma(u))*dx
L= dot(w,g_T)ds(2)+dot(w,g_T)dx

left_end_displacement = Constant((0.0, 0.0, 0.0))
bc = [DirichletBC(V, left_end_displacement, left)]
solve(a == L, u, bc)
plot(u, mode="displacement",interactive=True)

---

Hello ggrekas, Thanks for your help and time.

I have made the change suggested by you. However, the problem remains. I still get the same error for the following code:

from dolfin import *
d0 = 8.
d1 = 1.
d2 = 1.
n0 = int(d0*6)
n1 = int(d1*6)
n2 = int(d2*6)
o = Point(0, 0, 0)
l = Point(8, 1, 1)
Mesh = BoxMesh(o, l, n0, n1, n2)
Mesh.init()

V = VectorFunctionSpace(Mesh, "CG", 1)

tol = 1E-14 # tolerance

class Left(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 0.0)

class Right(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 8.0)

class Top(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[2], 1.0)

left = Left()
top = Top()
right = Right()

boundaries = FacetFunction("uint", Mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)

u = TrialFunction(V)
u = Function(V)           #EDIT
f = Constant((0.0, -1.9, 0.0))

E = 200*(10**9)
nu = 0.24

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

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

def sigma(v):
  return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().topological_dimension() )

T = TensorFunctionSpace(Mesh, 'CG', 1)
stress = project(sigma(u) ,T)      #EDIT
File("stress.pvd") << stress
M = stress.vector().array()
print M

ds = Measure("ds")[boundaries]

g_T= Constant((0.0, -8500000.0, 0.0))

w = TestFunction(V)
a = inner(grad(w),sigma(u))*dx
L = dot(w,g_T)*ds(2)+dot(w,g_T)*dx

left_end_displacement = Constant((0.0, 0.0, 0.0))
bc = [DirichletBC(V, left_end_displacement, left)] 
solve(a == L, u, bc)
plot(u, mode="displacement",interactive=True)

---

Hello again.

Before the definition of the bilinear and linear form, i.e. a and L
define u as
u = TrialFunction(V)

now before the calling of solve function u must be redefined as
u = Function(V)

One last change is on the stress computation because its definition
has been changed so
stress = project(sigma,T) must change to project(sigma(u),T)

I

---

Hello ggrekas, Thanks for your input.

I incorporated your changes (edited the code in my last comment above), however, now I get this error:

page@Jarvis:~$ python tensortrial.py
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving linear system of size 21609 x 21609 (PETSc Krylov solver).
[ 0.  0.  0. ...,  0.  0.  0.]
Traceback (most recent call last):
  File "tensortrial.py", line 69, in <module>
    solve(a == L, u, bc)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 297, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 321, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 82, in __init__
    cpp.LinearVariationalProblem.__init__(self, a, L, u, bcs)
  File "/usr/lib/python2.7/dist-packages/dolfin/cpp/fem.py", line 2412, in __init__
    _fem.LinearVariationalProblem_swiginit(self,_fem.new_LinearVariationalProblem(a, L, u, bcs))
RuntimeError: