Eigen solver syntax in dolfin-x

5

Hello sir,
I am attaching the dolfin code for the same problem as you have asked for it. I have to replicate the same using dolfin-x.

from dolfin import *
import numpy as np
import json
import matplotlib.pyplot as plt
import scipy as sp
%matplotlib inline
mesh = BoxMesh(Point(0.,0.,0.),Point(5,0.6,0.3), 25, 3, 2)
E, nu, rho = (2e11), (0.3), (7850)
mu,lmbda = E/2./(1+nu), E*nu/(1+nu)/(1-2*nu)
def sigma(v):
    dim = v.geometric_dimension()
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)
def eps(v):
    return sym(grad(v))
V = VectorFunctionSpace(mesh, 'CG', degree=2)
u = TrialFunction(V)
v = TestFunction(V)
left = CompiledSubDomain("near(x[0], 0) && on_boundary")
bc = DirichletBC(V, Constant((0.,0.,0.)), left)
k_form = inner(sigma(u),eps(v))*dx
m_form = rho*dot(u,v)*dx
K = PETScMatrix()
assemble(k_form, tensor=K)
M = PETScMatrix()
assemble(m_form, tensor=M)
bc.zero(M), bc.apply(K)
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.

N_eig = 40   # number of eigenvalues
print("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(N_eig)


L, B, H = 5., 0.3, 0.6 #Dimensions of the beam
# Exact solution computation
from scipy.optimize import root
from math import cos, cosh
falpha = lambda x: cos(x)*cosh(x)+1
alpha = lambda n: root(falpha, (2*n+1)*pi/2.)['x'][0]
file_results = XDMFFile("modal_analysis.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

# Extraction
eig_v = []
eig_l = []
for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)

    # 3D eigenfrequency
    freq_3D = sqrt(r)/2/pi

    # Beam eigenfrequency
    if i % 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(i/2)**2*sqrt(E*I_bend/(rho*B*H*L**4))/2/pi

    print("Solid FE: {0:8.5f} [Hz]   Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
    
    
    # Initialize function and assign eigenvector
    eigenmode = Function(V,name="Eigenvector "+str(i))
    eigenmode.vector()[:] = rx
    file_results.write(eigenmode, 0)
    eig_l.append(freq_3D)
    eig_v.append(rx)

The results are:

Computing 40 first eigenvalues...
Solid FE:  9.83051 [Hz]   Beam theory:  9.78457 [Hz]
Solid FE: 19.45146 [Hz]   Beam theory: 19.56914 [Hz]
Solid FE: 60.61639 [Hz]   Beam theory: 61.31884 [Hz]
Solid FE: 114.64296 [Hz]   Beam theory: 122.63769 [Hz]
Solid FE: 118.82391 [Hz]   Beam theory: 171.69454 [Hz]
Solid FE: 165.68699 [Hz]   Beam theory: 343.38908 [Hz]