Hello dear all, I am now modeling a 3D problem using the phase-field method. How to calculate eigenvalues and eigenvectors of the strain tensor in 3D? Can anyone guide me?
Best,
zgh,
You should post the problem you are working on or at least a minimum working example of what you have tried so others can help you. See here.
I can’t directly answer your question but I can provide a temporary work around.
Personally I encountered this same exact question a few months ago when I wanted to obtain the principle components of a symmetric stress tensor after essentially solving the Fenics tutorial for elasticity example.
-\nabla * \sigma(u) = f in \Omega
The closest I got was iterating through each of the vertices, getting the correct eigenvalues and eigenvectors using Numpy, but was unable to correctly project back onto the solution space. So as for doing this programmatically I am unable to help.
However, if you need this soon, the way I got around this was to export the projected solution \sigma(u) onto a TensorFunctionSpace and output as a “.xdmf”. I then opened this in Paraview and exported the tensor from the spreadsheet view as a “.csv”. Using another software (Matlab), I then could easily grab the stress tensor - calculate the components - and then generate a spreadsheet where I project the eigen values and vectors back onto the original solution space (to be view in ParaView). I found this to be optimal especially since my application involved 500k+ elements and even with MPI this became cumbersome to try in Python. This is my personal experience take it or leave it - I’m welcome to provide further detail if this is of interest.
-Sven
Dear Sven,
Thank you so much for your nice descriptions. More precisely in one part of the implementation of the model needs to eigenvalues and eigenvectors of strain tensor in case 3D. I found it only in case of 2D as follows:
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)
Unfortunately, I do not well know your way and I am looking for a direct way to solve this problem.
Hi Sven, im doing what you say and i can export my Stress Tensor as a csv file and calculate the principal Stress values in python and create a csv file whit these, but i dont know how project this into the original mesh in results, can you help me…
Do you plan to use the principal stresses only for visualization in ParaView or do you actually want to project the values back to a DG space in dolfin/dolfinx? I can advise on ParaView.
-Sven
Just for visualize it in paraview
Hope this helps (ParaView 5.7.0):
I exported the principal components and point coordinates (x,y,z) as a .csv file. I also generated (in legacy fenics) a mesh file of the domain using XDMFFile.
- Load both into ParaView.
- With the .csv selected, go to filters → Table To Points
- In the TableToPoints change the X, Y, and Z columns to the corresponding csv columns
- With the TableToPoints item selected go to filters → Point Dataset Interpolator
- Select the .csv as input, the mesh as source.
The corresponding object should be a solid color surface of your data. You can then select the coloring (options should be your other data columns in the csv) as normal.
-Sven
how can i generate the mesh file as XDMF??
ncarmona,
You should be able to follow the function description for XDMFFile. Example below is from that reference. This is what I used at the time with dolfin 2019.2.0.dev0.
with XDMFFile(mpi_comm_world(), 'name.xdmf') as xdmf:
xdmf.write(mesh)
-Sven
