# I meet a problem when running a linear variational problem and need to do a spectral decomposition of strain tensor **URL:** \ **Category:** Linear Algebra\ **Created:** [May 11, 2020, 3:39pm UTC](https://fenicsproject.discourse.group/t/i-meet-a-problem-when-running-a-linear-variational-problem-and-need-to-do-a-spectral-decomposition-of-strain-tensor/3169 "2020-05-11T15:39:56Z")\ **Posts on this page:** 6\ **Page:** 1
**Author:** ![Daniel](https://avatars.discourse-cdn.com/v4/letter/d/3ec8ea/32.png) [@Daniel](https://fenicsproject.discourse.group/u/Daniel)\ **Post date:** [May 11, 2020, 3:39pm UTC](https://fenicsproject.discourse.group/t/i-meet-a-problem-when-running-a-linear-variational-problem-and-need-to-do-a-spectral-decomposition-of-strain-tensor/3169/1 "2020-05-11T15:39:56Z")
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 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 "", line 1, in 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) ``` ---
**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\ **Post date:** [May 12, 2020, 7:10am UTC](https://fenicsproject.discourse.group/t/i-meet-a-problem-when-running-a-linear-variational-problem-and-need-to-do-a-spectral-decomposition-of-strain-tensor/3169/2 "2020-05-12T07:10:27Z")
Hi, you can have analytic expressions of eigenvalues in 2D directly in UFL instead of going through Sympy (I never tried replacing Sympy expressions with UFL but probably your problem comes from this step). If you then want to go to 3D, for which the above analytic expressions are not available, we have [a phase-field example](https://thelfer.github.io/mgis/web/mgis_fenics_phase_field.html) using MFront constitutive law code generator and FEniCS. ---
**Author:** ![Daniel](https://avatars.discourse-cdn.com/v4/letter/d/3ec8ea/32.png) [@Daniel](https://fenicsproject.discourse.group/u/Daniel)\ **Post date:** [May 13, 2020, 12:29pm UTC](https://fenicsproject.discourse.group/t/i-meet-a-problem-when-running-a-linear-variational-problem-and-need-to-do-a-spectral-decomposition-of-strain-tensor/3169/3 "2020-05-13T12:29:33Z")
many many thanks ,bleyerj. your mfont examples should be very help and i will carefully study it. some examples in your book “Numerical tours of continuum mechanics using FEniCS” is also very helpful for me as a beginner in FEniCS. thanks! ---
**Author:** ![Daniel](https://avatars.discourse-cdn.com/v4/letter/d/3ec8ea/32.png) [@Daniel](https://fenicsproject.discourse.group/u/Daniel)\ **Post date:** [May 14, 2020, 8:01am UTC](https://fenicsproject.discourse.group/t/i-meet-a-problem-when-running-a-linear-variational-problem-and-need-to-do-a-spectral-decomposition-of-strain-tensor/3169/4 "2020-05-14T08:01:43Z")
Hi, bleyerj. i opened the code of phase-field brittle fracture and find pycharm can not identify mgis.fenics module. I have done the following steps to install “MFrontGenericInterfaceSupport” 1 unzip field in `C:/MFrontGenericInterfaceSupport` (field from [this website](https://github.com/thelfer/MFrontGenericInterfaceSupport)) 2 type `cd /mnt/c/MFrontGenericInterfaceSupport` on terminal windows 3 type `cmake .` 4 type `cmake` 5 type `make` 6 type `make check` 7 type `make install` but after completing above steps, pycharm can not still identify mgis.fenics. i don’t know which step i wrongly make or lost. I have found in the “INSTALL-cmake.md” field, there are 2 options could be related but i do not know in which step i should use them. - `enable-python-bindings` : compiles the `Python` bindings. This requires the `Boost/Python` library to be available (default=OFF) - `enable-fenics-bindings` : compiles the `FEniCS` bindings. Those bindings are experimental and very limited. To use `MGIS` with `FEniCS` , usage of the `Python` bindings are encouraged. My OS is win10 and install ubuntu on it. ---
**Author:** ![xiaodong\_zhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/xiaodong_zhang/32/1135_2.png) [@xiaodong\_zhang](https://fenicsproject.discourse.group/u/xiaodong_zhang)\ **Post date:** [October 4, 2021, 8:44am UTC](https://fenicsproject.discourse.group/t/i-meet-a-problem-when-running-a-linear-variational-problem-and-need-to-do-a-spectral-decomposition-of-strain-tensor/3169/5 "2021-10-04T08:44:31Z")
Hi, Daniel. I don’t understand your procedure to get the positive and negative elastic energy, could you please explain it? ---
**Author:** ![dohertyw1](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dohertyw1/32/2561_2.png) [@dohertyw1](https://fenicsproject.discourse.group/u/dohertyw1)\ **Post date:** [July 26, 2022, 3:28pm UTC](https://fenicsproject.discourse.group/t/i-meet-a-problem-when-running-a-linear-variational-problem-and-need-to-do-a-spectral-decomposition-of-strain-tensor/3169/6 "2022-07-26T15:28:25Z")
Hi all, I have the same bug in the piece of code that Daniel posted. I need to calculate the eigenvalues of a 3x3 tensor and some link between sympy and ufl seems to be the best way to do it. Anyone have any idea how to help? Thanks, Will.