Is it possible to perform a spectral decomposition on the strain tensor with current ufl operators?

Question

I want to perform spectral decomposition on the strain tensor. Does FEniCS currently support such an operator? If not is it possible to write one?

Comments

Hi, do you intend to use this as part of defining a bilinear form? The following threads are relevant in any case, bad news, better news.

---

Hi MiroK,
Yes. I would like to separate the positive and negative principal components of the strain tensor after the spectral decomposition. In my variational form I wanted to penalize only the positive part of the strain tensor.
I had already looked at the better news and also this. But these cannot be used if I want to do a decomposition of a 2nd order strain tensor. I was hoping to write (even if it was challenging) a Python/C++ code to create an operator for the same. But from bad news it seems it is impossible. Please correct me if I am wrong.

Thanks
Vivek

---

Can you post relevant parts of the form?

---

# Constitutive relation
def stress( gradu, g):
    k0 = lmda + 2.*mu
    eps = sym( gradu )
    tr_pos_eps = conditional( gt(tr(eps), 0.), tr(eps) ,0) # tr(eps) == tr(eps_volumetric) as tr(eps_deviatoric) = 0
    tr_neg_eps = tr(eps) - tr_pos_eps
    dev_eps = dev(eps)
    return  k0*tr_neg_eps*Identity(space_dim) + g*(k0*Identity(space_dim)*tr_pos_eps + mu*dev_eps )

     # Stored Energy
def stored_energy( gradu, g ):
    k0 = lmda + 2.*mu
    eps = sym( gradu )
    tr_pos_eps = conditional( gt(tr(eps), 0.), tr(eps) ,0) # Shouldn't we divide this by space_dim?
    tr_neg_eps = tr(eps) - tr_pos_eps
    dev_eps = dev(eps) 
    return k0*pow(tr_neg_eps,2)/2. + g*(k0*pow(tr_pos_eps,2)/2. +  mu*inner(dev_eps,dev_eps) )

     # The penalization function
def penalize( phase_field ):
    return ( pow(1. - phase_field,2)  + 1.e-9 )

    # Set up the residuals
    R_u =  inner(stress( grad(uh) , penalize( dh ) ), grad( delta_u) )*dx + inner(source_u, delta_u)*dx 

Code Snippet could be found above. So instead of separating the strains into volumetric and deviatoric parts, I would like to do an eigenvalue decomposition.

Answer 1

Hi, see if the following helps. The idea here is to let sympy do the heavy lifting of coming up with the expression for eigenvalues and eigenvectors and then turning that to UFL. At the moment the conversion is via string and subsequent eval. You might want to build a better parser, especially because in 3d the expressions from sympy involve complex unit etc which do not have equivalent in UFL namespace.

from dolfin import *
import sympy as sp

# Compute the symbolic expression for eigenvalues by sympy
T = sp.Matrix(2, 2, lambda i, j: sp.Symbol('T[%d, %d]' % (i, j), real=True))
eig_expr = T.eigenvects()   # ((v, multiplicity, [w])...)

eigv = [e[0] for e in eig_expr]
eigw = [e[-1][0] for e in eig_expr]

eigv_expr = map(str, eigv)
eigw_expr = [[str(e[0]), str(e[1])] for e in eigw]

# UFL operator for eigenvalues of 2x2 matrix, a pair of scalars
def eigv(T): return map(eval, eigv_expr)

# UFL operator for eigenvectors of 2x2 matrix, a pair of vectors
def eigw(T): return [as_vector(map(eval, vec)) for vec in eigw_expr]

n = 16
mesh = UnitSquareMesh(n, n)

V = TensorFunctionSpace(mesh, 'CG', 1)
x = SpatialCoordinate(mesh)

v = as_vector((x[0]*x[1]+x[0], -x[0]*x[1]+x[1]))

t = sym(grad(v))

t = project(t, V)

v0, v1 = eigv(t)
w0, w1 = eigw(t)

# Eigenvalues okay
print sqrt(assemble(inner(det(t)-v0*v1, det(t)-v0*v1)*dx))
print sqrt(assemble(inner(tr(t)-(v0+v1), tr(t)-(v0+v1))*dx))

# Orthogonality of eigenvectors 
print sqrt(assemble(inner(Constant(0) - dot(w0, w1), Constant(0) - dot(w0, w1))*dx))

# Decomposition
e = dot(t, w0)
e0 = v0*w0
print sqrt(assemble(inner(e-e0, e-e0)*dx))

e = dot(t, w1)
e0 = v1*w1
print sqrt(assemble(inner(e-e0, e-e0)*dx))

Comments

Thanks a lot. I will try implementing the method.

---

Sure, and if possible please share your solution.

---

Hi Miroslav,

I tried implementing your suggestion of using sympy to obtain eigenvectors and eigenvalues, but I am running into troubles while reconstructing the strain tensor. I am getting an Arity error. In the following code I have attempted to simply reconstruct the strain using eigenvectors and eigenvalues but it throws an error. Code snippet:

from dolfin import *
import sympy as sp
import numpy as np
import numpy.linalg as LA

#----------------------------------------------------------------------#
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = 'uflacs'  
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# Compute the symbolic expression for eigenvalues by sympy
T = sp.Matrix(2, 2, lambda i, j: sp.Symbol('T[%d, %d]' % (i, j), real=True))
eig_expr = T.eigenvects()   # ((v, multiplicity, [w])...)

eigv = [e[0] for e in eig_expr]
eigw = [e[-1][0] for e in eig_expr]

eigv_expr = map(str, eigv)
eigw_expr = [[str(e[0]), str(e[1])] for e in eigw]

#----------------------------------------------------------------------#


#----------------------------------------------------------------------#
# Create UFL operators for the eigenvalues and eigenvectors

# UFL operator for eigenvalues of 2x2 matrix, a pair of scalars
def eigv(T): return map(eval, eigv_expr)

# UFL operator for eigenvectors of 2x2 matrix, a pair of vectors
def eigw(T): return [as_vector(map(eval, vec)) for vec in eigw_expr]

#----------------------------------------------------------------------#

spd = 2 # Spatial dimension

n = 16 # Size of mesh
mesh = UnitSquareMesh(n, n)

V = VectorFunctionSpace(mesh, 'CG', 1) # V = TensorFunctionSpace(mesh, 'CG', 1) -> Original
V_t = TensorFunctionSpace(mesh, 'CG', 1)

x = SpatialCoordinate(mesh)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0., -1.0))
T = Constant((0., 0.))

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)

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

lmbda = 10.0
mu = 1.0

def sigma(k):
    #return lmbda*(tr(eps(u)))*Identity(spd) +2*mu*eps(u)
    # To compare the answers
    t_e = eps(k)

    v0 , v1 = eigv(t_e)
    w0 , w1 = eigw(t_e)
    # v0, v1 are 'Sum objects
    # w0, w1 are 'ListTensor' objects

    # Creating tensor objects for eigenvalues
    vp = ([v0,0.0],[0.0,v1])

    vp = as_tensor(vp)

    #vp = project(vp,V_t) -> Gives Arity error

    # Creating tensor objects for eigenvectors
    wp = ([w0[0],w1[0]],[w0[1],w1[1]])
    wp = as_tensor(wp)
    # Taking the transpose of the eigenvector matrix
    wp_T = wp.T
    # Reconstructing the strain tensor
    strn = wp*vp*wp.T
    strn =sym(strn)

    # Using strn gives the following error:
    # ArityMismatch("Applying nonlinear operator {0} to expression depending on form argument {1}.".format(o._ufl_class_.__name__, t))
    # ArityMismatch: Applying nonlinear operator Sqrt to expression depending on form argument v_1.

    return lmbda*(tr(strn))*Identity(spd) +2.0*mu*strn

a = inner(sigma(u),eps(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds


# Compute solution
u = Function(V)
solve(a == L, u, bc)


V = FunctionSpace(mesh, 'CG', 1)

# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
print('min/max u:',
      u_magnitude.vector().array().min(),
      u_magnitude.vector().array().max())

I have already looked at your suggestion here but could not utilize it here. Any suggestions?

---

Hi Miroslav,

I attempted the following changes to test the reconstruction of the strain tensor and it seems to work in this code:

"""
The code performs an eignevalue decomposition of the strain tensor and
then reconstructs the strain tensor to compare the solutions.
"""

from dolfin import *
import sympy as sp
import numpy as np
import numpy.linalg as LA

#----------------------------------------------------------------------#
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = 'uflacs'  
#----------------------------------------------------------------------#


# Compute the symbolic expression for eigenvalues by sympy
T = sp.Matrix(2, 2, lambda i, j: sp.Symbol('T[%d, %d]' % (i, j), real=True))
eig_expr = T.eigenvects()   # ((v, multiplicity, [w])...)

eigv = [e[0] for e in eig_expr]
eigw = [e[-1][0] for e in eig_expr]

eigv = sp.Matrix(eigv) # Column vector
eigw = sp.Matrix(eigw) # Row has the elements
eigw = sp.Matrix([[eigw[0],eigw[2]],[eigw[1],eigw[3]]])
eigv = sp.Matrix([[eigv[0], 0.0],[0.0, eigv[1]]])

eig_dec = eigw*eigv*eigw.inv()

#print(T)
#print("--------------")
#print("eigw")
#print(eigw)
#print("--------------")
#print("eigv")
#print(eigv)
#print("--------------")
#print("eigw.inv()")
#print(eigw.inv())
#print("--------------")
#print("--------------")
#print("eig_dec")
#print(sp.simplify(eig_dec))


eigv_expr = map(str, eigv)
eigw_expr = map(str, eigw)


# UFL operator for eigenvalues of 2x2 matrix, a pair of scalars
def eigv(T): return map(eval, eigv_expr)

# UFL operator for eigenvectors of 2x2 matrix, a pair of vectors
def eigw(T): return map(eval, eigw_expr)


n = 16
mesh = UnitSquareMesh(n, n)

V = TensorFunctionSpace(mesh, 'CG', 1)
x = SpatialCoordinate(mesh)

v = as_vector((x[0]*x[1]+x[0], -x[0]*x[1]+x[1]))

t = sym(grad(v))

t = project(t, V)

v00, v01, v10, v11 = eigv(t)

w00, w01, w10, w11 = eigw(t)

w = ([w00,w01],[w10,w11])
w = as_tensor(w)

v = ([v00,v01],[v10,v11])
v = as_tensor(v)

strn = w*v*inv(w)

# Comparing the strain obtained using eigenvalue decomposition with the original strain.
print sqrt(assemble(inner(t-strn, t-strn)*dx))

print sqrt(assemble(inner(det(t)-v[0,0]*v[1,1], det(t)-v[0,0]*v[1,1])*dx))

print sqrt(assemble(inner(tr(t)-tr(strn), tr(t)-tr(strn))*dx))

#Looks ok

I then tried the same implementation with a 2-D elasticity problem and again ran into the Arity issue.

from dolfin import *
import sympy as sp
import numpy as np
import numpy.linalg as LA

#----------------------------------------------------------------------#
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = 'uflacs'  
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# Compute the symbolic expression for eigenvalues by sympy
T = sp.Matrix(2, 2, lambda i, j: sp.Symbol('T[%d, %d]' % (i, j), symmetric = True, real=True))
eig_expr = T.eigenvects()   # ((v, multiplicity, [w])...)

eigv = [e[0] for e in eig_expr]
eigw = [e[-1][0] for e in eig_expr]

eigv = sp.Matrix(eigv) # Column vector
eigw = sp.Matrix(eigw) # Row has the elements
eigw = sp.Matrix([[eigw[0],eigw[2]],[eigw[1],eigw[3]]])
eigv = sp.Matrix([[eigv[0], 0.0],[0.0, eigv[1]]])

eigv_expr = map(str, eigv)
eigw_expr = map(str, eigw)

#----------------------------------------------------------------------#


#----------------------------------------------------------------------#
# Create UFL operators for the eigenvalues and eigenvectors

# UFL operator for eigenvalues of 2x2 matrix, a pair of scalars
def eigv(T): return map(eval, eigv_expr)

# UFL operator for eigenvectors of 2x2 matrix, a pair of vectors
def eigw(T): return map(eval, eigw_expr)

#----------------------------------------------------------------------#

spd = 2 # Spatial dimension

n = 16 # Size of mesh
mesh = UnitSquareMesh(n, n)

V = VectorFunctionSpace(mesh, 'CG', 1) # V = TensorFunctionSpace(mesh, 'CG', 1) -> Original
V_t = TensorFunctionSpace(mesh, 'CG', 1)

x = SpatialCoordinate(mesh)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0., -1.0))
T = Constant((0., 0.))

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)

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


def strn(u):
    t_e = eps(u)

    v00, v01, v10, v11 = eigv(t_e)

    w00, w01, w10, w11 = eigw(t_e)

    wp = ([w00,w01],[w10,w11])
    wp = as_tensor(wp)

    vp = ([v00,v01],[v10,v11])
    vp = as_tensor(vp)

    return wp*vp*inv(wp)



lmbda = 10.0
mu = 1.0

def sigma(u):
    #return lmbda*(tr(eps(u)))*Identity(spd) +2*mu*eps(u)
    # To compare the answers

    # Using strn instead of eps(u) gives the following error:
    # ArityMismatch("Applying nonlinear operator {0} to expression depending on form argument {1}.".format(o._ufl_class_.__name__, t))
    # ArityMismatch: Applying nonlinear operator Sqrt to expression depending on form argument v_1.
    # Here : https://fenicsproject.org/qa/6882/problem-with-mixed-space-formulation
    # It is mentioned that the object to be projected must not be defined in terms of trial/test functions
    # Seems irrelevant but the only one that is close enough

    return lmbda*(tr(strn(u)))*Identity(spd) +2.0*mu*strn(u)

a = inner(sigma(u),eps(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds


# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Not adding the following line gives an error
# Why do we need to add this here?
V = FunctionSpace(mesh, 'CG', 1)

# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
#plot(u_magnitude, 'Displacement magnitude')
# Any suggestions?
print('min/max u:',
      u_magnitude.vector().array().min(),
      u_magnitude.vector().array().max())

Error message

in nonlinear_operator
raise ArityMismatch("Applying nonlinear operator {0} to expression depending on form argument {1}.".format(o._ufl_class_.name, t))
ArityMismatch: Applying nonlinear operator Sqrt to expression depending on form argument v_1.

Any suggestions?

---

Basically function sigma above is nonlinear as a consequence of eigenpair computations. You need to set up the problem as nonlinear with u as a Function.

---

Thanks! It works.

---

Update:
Even though it worked, i.e. no Arity errors. My Newton solver is no longer able to solve the problem. The errors shown are : r(abs)=-nan; r(rel)=-nan. I am looking into the possible sources of errors.

---

I figured out the error. It was the incorrect initialization of the solution function.

Thanks a lot for the help MiroK

---

Perfect, could you update the code you posted in the comment such that it works? Thanks.

---

Added it as an answer. Thanks!

Answer 2

I am answering my own question so that a working code is available:

Code for eigenvalue decomposition:

from dolfin import *
import sympy as sp
import numpy as np
import numpy.linalg as LA

#----------------------------------------------------------------------#
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = 'uflacs'  
#----------------------------------------------------------------------#

#----------------------------------------------------------------------#
# Compute the symbolic expression for eigenvalues by sympy
T = sp.Matrix(2, 2, lambda i, j: sp.Symbol('T[%d, %d]' % (i, j), symmetric =True, real=True))
eig_expr = T.eigenvects()   # ((v, multiplicity, [w])...)
#eig_expr = T.LDLdecomposition()

eigv = [e[0] for e in eig_expr]
eigw = [e[-1][0] for e in eig_expr]

eigv = sp.Matrix(eigv) # Column vector
eigw = sp.Matrix(eigw) # Row has the elements
eigw = sp.Matrix([[eigw[0],eigw[2]],[eigw[1],eigw[3]]])
eigv = sp.Matrix([[eigv[0], 0.0],[0.0, eigv[1]]])

eigv_expr = map(str, eigv)
eigw_expr = map(str, eigw)

#----------------------------------------------------------------------#


#----------------------------------------------------------------------#
# Create UFL operators for the eigenvalues and eigenvectors

# UFL operator for eigenvalues of 2x2 matrix, a pair of scalars
def eigv(T): return map(eval, eigv_expr)

# UFL operator for eigenvectors of 2x2 matrix, a pair of vectors
def eigw(T): return map(eval, eigw_expr)

#----------------------------------------------------------------------#

spd = 2 # Spatial dimension

n = 16 # Size of mesh
mesh = UnitSquareMesh(n, n)

V = VectorFunctionSpace(mesh, 'CG', 1) # V = TensorFunctionSpace(mesh, 'CG', 1) -> Original

x = SpatialCoordinate(mesh)

# Define variational problem
du = TrialFunction(V)
delta_u = TestFunction(V)

# set up solution functions
uh = Function(V,name='displacement')
# The way the eigenvalues are computed we cannot allow a constant value of u at start

uh_array = uh.vector().array()
uh_array = np.random.rand(len(uh_array))
uh.vector()[:] = uh_array


force = Constant((0., -1.0))
T = Constant((0., 0.))

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0)), clamped_boundary)

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


def strn(u):
    t = sym(grad(u))

    v00, v01, v10, v11 = eigv(t)

    w00, w01, w10, w11 = eigw(t)

    w = ([w00,w01],[w10,w11])
    w = as_tensor(w)

    v = ([v00,v01],[v10,v11])

    v = as_tensor(v)


    return w*v*inv(w)
# Positive strain
def strn_p(u):
    t = sym(grad(u))
    v00, v01, v10, v11 = eigv(t)

    v00 = conditional(gt(v00,0.0),v00,0.0)
    v11 = conditional(gt(v11,0.0),v11,0.0)

    w00, w01, w10, w11 = eigw(t)
    wp = ([w00,w01],[w10,w11])
    wp = as_tensor(wp)
    vp = ([v00,v01],[v10,v11])
    vp = as_tensor(vp)  
    return wp*vp*inv(wp)
# Negative strain
def strn_n(u):
    t = sym(grad(u))
    v00, v01, v10, v11 = eigv(t)
    v00 = conditional(lt(v00,0.0),v00,0.0)
    v11 = conditional(lt(v11,0.0),v11,0.0)
    w00, w01, w10, w11 = eigw(t)
    wn = ([w00,w01],[w10,w11])
    wn = as_tensor(wn)
    vn = ([v00,v01],[v10,v11])
    vn = as_tensor(vn)  
    return wn*vn*inv(wn)

lmbda = 10.0
mu = 1.0

def sigma(u):

    #return lmbda*(tr(eps(u)))*Identity(spd) +2*mu*eps(u)

    # To compare the answers

    return lmbda*(tr(strn_p(u)+strn_n(u)))*Identity(spd) +2.0*mu*(strn_p(u)+strn_n(u))




F_u = inner(sigma(uh),grad(delta_u))*dx + dot(force, delta_u)*dx

J_u = derivative(F_u, uh, du )

# Compute solution
problem = NonlinearVariationalProblem(F_u, uh, bc, J_u)
solver  = NonlinearVariationalSolver(problem)

solver.solve()

V = FunctionSpace(mesh, 'CG', 1)

# Compute magnitude of displacement
u_magnitude = sqrt(dot(uh, uh))
u_magnitude = project(u_magnitude, V)
#plot(u_magnitude, 'Displacement magnitude')
# Any suggestions?
print('min/max u:',
      u_magnitude.vector().array().min(),
      u_magnitude.vector().array().max())

#----------------------------------------------------------------------#
# A simple 2-D elasticity problem ends here.
# The solution obtained was:
# ('min/max u:', -5.0591841289798924e-05, 0.9251112919755422)
#
# Reconstructing the strain gives:
#('min/max u:', -5.0591841289772903e-05, 0.92511129197543296)
# Obtaining the values manually and solving gives me the same answer.
#----------------------------------------------------------------------#

Comments

Cool, is there a need for inv in str_n and str_p? The eigenvectors should be orthonormal and therefore the inverse could be computed by transpose and in your case the expression for it defined directly. In any case, nicely done.

---

Hi Miroslav,

Using the transpose instead of inverse was giving erroneous results. I am not exactly sure why. One possibility could be that the eigenvalue decomposition does not lead to orthogonal vectors. If could find the reason and the solution, I would update the code. I am looking into it so that the computations could be faster.

Thanks