Eigen vectors in dolfinx

1

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,
                 nabla_grad, dx)

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):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)


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
Screenshot 2022-01-06 083818
Can you please help me in computing and writing eigen vectors in xdmf format

2

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
3

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
image