Sir, Can you rerun this code.
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,assemble_matrix, form)
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= form(inner(sigma(u), epsilon(v)) * dx)
K = assemble_matrix(k_form, bcs=[bc])
K.assemble()
m_form = form(rho*dot(u, v)*dx)
M = 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]
import dolfinx.io
xdmf = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "eigenvector.xdmf", "w")
xdmf.write_mesh(mesh)
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))
ur = dolfinx.fem.Function(V)
ur.vector.array[:] = vr
xdmf.write_function(ur.copy())
eig_v.append(vr.copy())
ef += 1
xdmf.close()
I am getting white contour of eigen vector with Paraview 5.10.0
