hello, I am now modelling a 2D-problem using the phase field method. I need to decompose strain tensor into positive (tension) parts and negatives parts. but when i run the following code, i get the error listed in the following.
Thanks in advance
Traceback (most recent call last):
File"/mnt/c/Users/Danie/OneDrive/PhaseFieldFEniCS/PhaseFieldFEniCS/PhaseField.py", line 169, in <module>
E_du = inner(sigma(u),grad(v))*dx
File"/mnt/c/Users/Danie/OneDrive/PhaseFieldFEniCS/PhaseFieldFEniCS/PhaseField.py", line 155, in sigma
return ((1-pold)**2)*(lmbda * tr_p * Identity(2) + 2.0 * mu * strn_p(u))+lmbda * tr_n * Identity(2)+ 2.0 * mu * strn_n(u)
File"/mnt/c/Users/Danie/OneDrive/PhaseFieldFEniCS/PhaseFieldFEniCS/PhaseField.py", line 124, in strn_p
v00, v01, v10, v11 = eigv(t)
File "<string>", line 1, in <module>
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 605, in sqrt
return _mathfunction(f, Sqrt)
File "/usr/lib/python3/dist-packages/ufl/operators.py", line 594, in _mathfunction
f = as_ufl(f)
File "/usr/lib/python3/dist-packages/ufl/constantvalue.py", line 472, in as_ufl
" to any UFL type." % str(expression))
ufl.log.UFLValueError: Invalid type conversion: T[0, 0]**2 - 2*T[0, 0]*T[1, 1] + 4*T[0, 1]*T[1, 0] + T[1, 1]**2 can not be converted to any UFL type.
the following is code.
import dx
import sympy
import numpy
from dolfin import *
#----------------------------------------------------------------------#
# 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 = sympy.Matrix(2, 2, lambda i, j: sympy.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 = sympy.Matrix(eigv) # Column vector
eigw = sympy.Matrix(eigw) # Row has the elements
eigw = sympy.Matrix([[eigw[0],eigw[2]],[eigw[1],eigw[3]]])
eigv = sympy.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)
#----------------------------------------------------------------------#
#----------------------------------------------------------------------#
# Import mesh from .xml
mesh = Mesh('mesh.xml')
# Introduce manually the material parameters
Gc = 2.7
l = 0.015
lmbda = 121.1538e3
mu = 80.7692e3
#----------------------------------------------------------------------#
#----------------------------------------------------------------------#
# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
p, q = TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)
#----------------------------------------------------------------------#
#----------------------------------------------------------------------#
# set up solution functions
u,unew, uold = Function(W,name='displacement'), Function(W,name='displacement'), Function(W,name='displacement')
pnew, pold, Hold = Function(V,name='phi'), Function(V,name='phi'), Function(V,name='phi')
# The way the eigenvalues are computed we cannot allow a constant value of u at start
unew_array = numpy.array(unew.vector())
unew_array = numpy.random.rand(len(unew_array))
unew.vector()[:] = unew_array
uold_array = numpy.array(uold.vector())
uold_array = numpy.random.rand(len(uold_array))
uold.vector()[:] = uold_array
u_array = numpy.array(u.vector())
u_array = numpy.random.rand(len(u_array))
u.vector()[:] = u_array
#----------------------------------------------------------------------#
#----------------------------------------------------------------------#
# Boundary conditions and control
top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")
load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(0), load, top)
bc_u = [bcbot, bctop]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
#----------------------------------------------------------------------#
#----------------------------------------------------------------------#
#performing strain tensor spectral decomposition
def epsilon(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)
#----------------------------------------------------------------------#
# Define stress Constituive functions
def sigma(u):
#return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))
# To compare the answers
tr_p=0.5 * (tr(epsilon(u)) + abs(tr(epsilon(u))))
tr_n=0.5 * (tr(epsilon(u)) - abs(tr(epsilon(u))))
return ((1-pold)**2)*(lmbda * tr_p * Identity(2) + 2.0 * mu * strn_p(u))+lmbda * tr_n * Identity(2)+ 2.0 * mu * strn_n(u)
#----------------------------------------------------------------------#
#----------------------------------------------------------------------#
# define positive elastic energy density and introduce a strain-history field
def psi(u):
return 0.5*(lmbda+mu)*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2+\
mu*inner(dev(epsilon(u)),dev(epsilon(u)))
def H(uold,unew,Hold):
return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)
#----------------------------------------------------------------------#
#----------------------------------------------------------------------#
#Writing governing equations in weak form
E_du = inner(sigma(u),grad(v))*dx
E_phi = (Gc*l*inner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))\
*inner(p,q)-2.0*H(uold,unew,Hold)*q)*dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi), rhs(E_phi), pnew)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)