Information about SLEPcEigenSolver

Linear Algebra
1

Hello,

I am trying to compute the principal strains and their associated principal direction in the case of linear elasticity. I noticed the SLEPcEigenSolver which can compute the eigenvalues and eigenvectors of a given matrix.

Could it be used for this purpose?

Here is a MWE for the usage of the solver in the case of a 2D rectangle under traction, however I think I am computing the eigenvalues and eigenvectors of the stiffness matrix.
Is there a way to pass the strain tensor instead?

Thanks in advance

from dolfin import *
import mshr as mr

# --- Mesh
L, H = 1, 0.1
nel = 75
geom = mr.Rectangle(Point(0, 0), Point(L, H))
mesh = mr.generate_mesh(geom, nel)
ndim = mesh.topology().dim()

# --- Material Parameter
E = Constant(70e3)
nu = Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
k = lmbda + mu*2/3.

# --- Function Spaces
Vs = FunctionSpace(mesh, 'DG', 0)
V = VectorElement("Lagrange", mesh.ufl_cell(), 1)
V_u = FunctionSpace(mesh, V)

u = Function(V_u)
du = TrialFunction(V_u)
v = TestFunction(V_u)
u.rename('u', 'u')

S3 = TensorFunctionSpace(mesh, 'DG', 0)
Eps = Function(S3)
Eps.rename('E', 'E')

# ---Boundaries
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
left = Left()

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

def NodeL(x):
    return near(x[0], 0) and near(x[1], 0.)

# ---Measures
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 3)

# ---Displacement BC
#u_R = Expression("t*(0.5-x[1]/H)", H=H, t=0.05, degree=1)
u_R = Expression("t", t=0.05, degree=0)
bcu_l = DirichletBC(V_u.sub(0), Constant(0.), boundaries, 1)
bcu_r = DirichletBC(V_u.sub(0), u_R, boundaries, 3)
bcu_L = DirichletBC(V_u.sub(1), Constant(0.), NodeL, method = 'pointwise')
bcu = [bcu_r, bcu_l, bcu_L]

# ---Definitions
def eps(u):
    return sym(grad(u))

def sigma(u):
    return 2. * mu * (eps(u)) + lmbda * tr(eps(u)) * Identity(2)

# --- Eigenproblem
k_form = inner(sigma(v), eps(du)) *dx
K = PETScMatrix()
assemble_system(k_form, l_form, bcu, A_tensor=K)
eigensolver = SLEPcEigenSolver(K)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 1e-12

N_eig = 2
eigensolver.solve(N_eig)

eI_, c, nI_, cx = eigensolver.get_eigenpair(0)
eI = Function(Vs, name = 'evI')
eI.vector()[:] = eI_
nI = Function(V_u, name = "nI")
nI.vector()[:] = nI_

eII_, c, nII_, cx = eigensolver.get_eigenpair(1)
eII = Function(Vs, name = 'evII')
eII.vector()[:] = eII_
nII = Function(V_u, name = "nII")
nII.vector()[:] = nII_

File('Linear_elasticity/' + 'evI.pvd') << eI
File('Linear_elasticity/' + 'evII.pvd') << eII
File('Linear_elasticity/' + 'nI.pvd') << nI
File('Linear_elasticity/' + 'nII.pvd') << nII
2

Hi, in principle you could but this really not what SLEPc is intended for. It is dedicated to solving “global” eigenvalue problems like modal or buckling analysis of a structure.
In 2D you have explicit analytical expressions of eigenvalues and eigenvectors that you can implement symbolically using UFL.
Otherwise, if your use case is only for postprocessing purposes, I suggest looping over the strain values when projected on a suitable function space and compute analytically or numerically (using Numpy or Scipy eig function) the eigenvalues/eigenvectors. You could also use Paraview’s filter to compute eigenvalues.

2 Likes
3

Thank you very much for the reply.

I was trying alternative routes since I foud that the UFL implementation for the spectral decomposition of the strain tensor presents some issues when damage (i.e.: phase-field method for brittle fracture) affects only the positive part of the strain energy.

Even when using the Numpy/SciPy eig function to eavluate the eigenvalues and eigenvectors, then projecting them over DG scalar/vector functions, I found it to be affected by similar problems as well.

I was looking for some other possible implementation methods to recover the eigenvectors since they assume some erratic behaviour.

For context, here are the UFL and the numerical implementations for the evaluation of both the eigenvalues and eigenvectors I am currently using.

def eigvalues(T):
    TI, TII = (tr(T) + sqrt(tr(T)**2 - 4*det(T)))/2, (tr(T) - sqrt(tr(T)**2 - 4*det(T)))/2
    return conditional(TI >= TII, as_vector([TI, TII]), as_vector([TII, TI]))

def eigvectors(T):
    tol = DOLFIN_EPS
    TI, TII = eigvalues(T)[0], eigvalues(T)[1]
    Txx, Txy, Tyy = T[0,0], T[0,1], T[1,1]
    theta_I = uf.atan_2(TI - Txx, Txy)
    theta_II = theta_I + np.pi/2
    vI = conditional(gt(abs(Txy), tol), as_vector([cos(theta_I), sin(theta_I)]), (conditional(ge(Txx, Tyy), as_vector([1., 0.]), as_vector([0., 1.]))))
    vII = conditional(gt(abs(Txy), tol), as_vector([cos(theta_II), sin(theta_II)]), (conditional(ge(Txx, Tyy), as_vector([0., 1.]), as_vector([1., 0.]))))
    return [vI, vII]

S3 = TensorFunctionSpace(mesh, 'DG', 0)
S1 = FunctionSpace(mesh, 'DG', 0)

def eigh_np2dlf(ep):
    ep_proj = project(ep, S3)
    evI = Function(S1)
    evII = Function(S1)
    nIx = Function(S1)
    nIy = Function(S1)
    nIIx = Function(S1)
    nIIy = Function(S1)
    for n, cell in enumerate(cells(mesh)):
        eps_loc = ep_proj.vector()[n*4:(n+1)*4]
        eps_array = np.array(eps_loc).reshape(2,2)
        ev_np, n_np = np.linalg.eigh(eps_array)
        if ev_np[0] >= ev_np[1]:
            evI.vector()[n] = ev_np[0]
            nIx.vector()[n] = n_np[:, 0][0]
            nIy.vector()[n] = n_np[:, 0][1]
            evII.vector()[n] = ev_np[1]
            nIIx.vector()[n] = n_np[:, 1][0]
            nIIy.vector()[n] = n_np[:, 1][1]
        else:
            evI.vector()[n] = ev_np[1]
            nIx.vector()[n] = n_np[:, 1][0]
            nIy.vector()[n] = n_np[:, 1][1]
            evII.vector()[n] = ev_np[0]
            nIIx.vector()[n] = n_np[:, 0][0]
            nIIy.vector()[n] = n_np[:, 0][1]

    nI = as_vector([nIx, nIy])
    nII = as_vector([nIIx, nIIy])
    return evI, evII, nI, nII
4

So what problem do you encounter exactly with such implementations ?

5

Sorry for the late (and probably too long) reply which drifts a bit from the initial post as it delves a little more into the problem I am trying to solve.
I tried my best to be as brief as i could while at the same time trying to provide all the useful informations which could help identify the problems in the evaluation of the eigenvalues/eigenvectors.

Here is the test setup I considered

image
image761Ă—154 1.21 KB

First, I started with the standard phase-field formulation for brittle materials (alternate minimization algorithm), only saving the eigenvalues and eigenvectors at the end of each iteration without using them into the formulation.
While in the pure linear elastic phase both method provide the same results, during the load increment where damage nucleates, the eigenvectors results from the two methods differs substantially (first componenet of the first eigenvector)

Eventually, I used the “No-tension material” phase-field formulation proposed by Freddi, Carfagni (Redirecting) where the strain energy is splitted into its positive and negative parts, and the damage is applied only to the positive one.
Here, even if the NumPy results were more accurate in the previous case, as the fracture starts to localize the eigenvectors changes and the problem does not reach convergence.

np_nIx
np_nIx1156Ă—382 58.1 KB

On the other hand, the UFL formulation reaches convergence even if the eigenvectors behave in a very strange way. The results of the first component of the first egienvector is reported below at it = 15 and at it = end of load increment.
UFL_nIx
UFL_nIx939Ă—476 53 KB

I am kinda confused by the discrepancy between the results as both methods utilize the same starting strain tensor from which the eigenvalues and eigenvectors are evaluated.
While I understand that in order to use NumPy I am performing a projection which could introduce some errors as well as that NumPy uses a different algorithm to numerically solve the eigenproblem, the differences between the two methods are quite significatives.

Again sorry for the super long reply, and if you reached this far thank you very much for your time reading through this.
(If i missed something, which is very likely, let me know and i will try to fix it)