# Eigen vectors in dolfinx

Hi, I am solving an eigen value analysis/modal analysis problem. I want to compute `eigen vectors` in dolfin-x. The MWE for the same has been attached. I believe, the `vr` PETSc vectors contain the eigenvector.

``````import sys
import dolfinx
import numpy as np
import slepc4py
from scipy.optimize import root
from mpi4py import MPI
from petsc4py import PETSc
from slepc4py import SLEPc

from dolfinx.mesh import CellType, locate_entities_boundary, create_box
from dolfinx.fem import (Constant, DirichletBC, Function, VectorFunctionSpace,
locate_dofs_topological)
from ufl import (TrialFunction, TestFunction, Identity, inner, dot, nabla_div,

slepc4py.init(sys.argv)

L, B, H = 5., 0.3, 0.6

mesh = create_box(
MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([L, H, B])],
[25, 3, 2],
cell_type=CellType.tetrahedron)

E, nu = (2E11), (0.3)
rho = (7850)
mu = E/2./(1+nu)
lambda_ = E*nu/(1+nu)/(1-2*nu)

def epsilon(u):

def sigma(u):
return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) \
+ 2*mu*epsilon(u)

V = VectorFunctionSpace(mesh, ("CG", 2))

u = TrialFunction(V)
v = TestFunction(V)

def clamped_boundary(x):
return np.isclose(x[0], 0)

fdim = mesh.topology.dim - 1
boundary_facets = locate_entities_boundary(mesh, fdim, clamped_boundary)

u_D = Function(V)
with u_D.vector.localForm() as loc:
loc.set(0)

bc = DirichletBC(u_D, locate_dofs_topological(V, fdim, boundary_facets))

T = Constant(mesh, (PETSc.ScalarType(0), PETSc.ScalarType(0),
PETSc.ScalarType(0)))
k_form = inner(sigma(u), epsilon(v))*dx
m_form = rho*dot(u, v)*dx

K = dolfinx.fem.assemble_matrix(k_form, bcs=[bc])
K.assemble()

M = dolfinx.fem.assemble_matrix(m_form, bcs=[bc])
M.assemble()

shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType('sinvert')  # spectral transform
shift.setShift(0.0)  # spectral shift
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setDimensions(nev=50)
eigensolver.setProblemType(2)
eigensolver.setST(shift)
eigensolver.setType('krylovschur')  # Solver type
eigensolver.setWhichEigenpairs(7)  # Target Real
eigensolver.setConvergenceTest(2)
eigensolver.setOperators(K, M)

eps_type = eigensolver.getType()
print("Solution method: %s" % eps_type)
problem_type = eigensolver.getProblemType()
print("Problem type: %s" % problem_type)
which_eigenpairs = eigensolver.getWhichEigenpairs()
print("Eigenpairs type: %s" % which_eigenpairs)

eigensolver.solve()
evs = eigensolver.getConverged()
vr, vi = K.createVecs()

# Exact solution computation
def falpha(x):
return np.cos(x) * np.cosh(x) + 1

def alpha(n):
return root(falpha, (2*n + 1)*np.pi/2)['x'][0]

eig_v=[]
ef = 0
print("Number of converged eigenpairs %d" % evs)
if evs > 0:
for i in range(evs):
e_val = eigensolver.getEigenpair(i, vr, vi)

if (~np.isclose(e_val.real, 1.0)):
# Calculation of eigenfrequency from real part of eigenvalue
freq_3D = np.sqrt(e_val.real)/2/np.pi

# Beam eigenfrequency
if ef % 2 == 0:
# exact solution should correspond to weak axis bending
I_bend = H*B**3/12.
else:
# exact solution should correspond to strong axis bending
I_bend = B*H**3/12.

freq_beam = alpha(ef/2)**2*np.sqrt(E*I_bend/(rho*B*H*L**4))/2/np.pi

print(
"Solid FE: {0:8.5f} [Hz] "
"Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
eig_v.append(vr)
ef += 1
print(np.array(eig_v).T)
``````

With this code getting this result

``````array([[ 7.87128161e-05,  7.87128161e-05,  7.87128161e-05, ...,
7.87128161e-05,  7.87128161e-05,  7.87128161e-05],
[-4.79422279e-05, -4.79422279e-05, -4.79422279e-05, ...,
-4.79422279e-05, -4.79422279e-05, -4.79422279e-05],
[ 2.09577501e-06,  2.09577501e-06,  2.09577501e-06, ...,
2.09577501e-06,  2.09577501e-06,  2.09577501e-06],
...,
[-3.39323655e-03, -3.39323655e-03, -3.39323655e-03, ...,
-3.39323655e-03, -3.39323655e-03, -3.39323655e-03],
[-1.81604351e-02, -1.81604351e-02, -1.81604351e-02, ...,
-1.81604351e-02, -1.81604351e-02, -1.81604351e-02],
[-2.21205070e-03, -2.21205070e-03, -2.21205070e-03, ...,
-2.21205070e-03, -2.21205070e-03, -2.21205070e-03]])
``````

but the correct values from dolfin code for same problem statement is

Can you please help me in computing and writing eigen vectors in `xdmf` format

You have a mistake in the line

``````eig_v.append(vr)
``````

Try changing it to

``````eig_v.append(vr.copy())
``````

Note that `dolfin` and `dolfinx` handle Dirichlet βrowsβ differently. `dolfinx` will place β1β on the diagonal for Dirichlet rows. In `dolfin` the value depends on the number of connected cells.

3 Likes

Thanks a lot sir for the correction. I tried with this. But results are not matching with dolfin. This is the result that I a getting for dolfinx