Hello all.
Fenicsx finds the slepc module at import time, but when
eigs = slepc.EPS().create(mesh.comm)
is encountered, the program crashes with the following error:
Traceback (most recent call last):
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection2.py", line 127, in <module>
eigs = slepc.EPS().create(mesh.comm)
AttributeError: module 'slepc4py' has no attribute 'EPS'
This is the code I used.
# Generate unit cell of helical structure (no periodic boundary in this case).
import sys
import numpy as np
from mpi4py import MPI as nMPI
from petsc4py import PETSc
import slepc4py as slepc
import meshio
import ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx import fem
from dolfinx import io
import gmsh
l = 6.0 # Cell length
r = 6.0 # Cell radius
rh = 2.0 # Helix radius
rw = 0.2 # wire radius
npts = 100 # number of points of spline curve
eps = 1.0e-5
Gamma = 0 + 6.283j
comm = nMPI.COMM_WORLD
mpiRank = comm.rank
modelRank = 0 # process where GMSG runs
if mpiRank == modelRank:
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 1)
gmsh.model.add("Helix Cell")
c = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, l, r, 1, 2*np.pi)
dt = 2 * np.pi / npts
p = []
for n in range(npts+3):
theta = 2 * np.pi * (n-1) / npts
kk = gmsh.model.occ.addCircle(rh * np.cos(theta), rh * np.sin(theta), (n -1) * l / npts, rw)
p.append(gmsh.model.occ.addCurveLoop([kk], tag=-1))
print(p)
s = gmsh.model.occ.addThruSections(p, tag=-1)
print(s)
v = gmsh.model.occ.cut([(3,c)], s, tag=100, removeObject=True, removeTool=True)
print("v={}".format(v))
for n in range(npts+3):
gmsh.model.occ.remove([(1,p[n])], recursive=True)
gmsh.model.occ.synchronize()
bb = gmsh.model.getEntities(dim=3)
pt = gmsh.model.getBoundary(bb, combined=True, oriented=False, recursive=False)
print(pt)
# Boundary conditions
Bot = []
Top = []
PEC = []
for bb in pt:
CoM = gmsh.model.occ.getCenterOfMass(bb[0], bb[1])
print(CoM)
if(np.abs(CoM[2]) < eps):
Bot.append(bb[1])
elif(np.abs(CoM[2]-l) < eps):
Top.append(bb[1])
else:
PEC.append(bb[1])
print(Bot, Top, PEC)
gmsh.model.mesh.setPeriodic(2, Top, Bot, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, l, 0, 0, 0, 1])
gmsh.model.addPhysicalGroup(3, [100], 1) # Helix volume
gmsh.model.setPhysicalName(3, 1, "HelixVolume")
gmsh.model.addPhysicalGroup(2, Top, 1) # Top BC (slave periodic BC)
gmsh.model.addPhysicalGroup(2, Bot, 2) # Bottom BC (master periodic BC)
gmsh.model.addPhysicalGroup(2, PEC, 3) # PEC boundary
gmsh.model.setPhysicalName(2, 1, "TopPBC")
gmsh.model.setPhysicalName(2, 2, "BotPBC")
gmsh.model.setPhysicalName(2, 3, "PEC")
gmsh.model.mesh.setSize([(0,108),(0,109)], 0.05)
gmsh.option.setNumber('Mesh.MeshSizeMin', 0.001)
gmsh.option.setNumber('Mesh.MeshSizeMax', 1.2)
gmsh.option.setNumber('Mesh.Algorithm', 6) #1=Mesh Adapt, 2=Auto, 3=Initial mesh only, 5=Delaunay, 6=Frontal-Delaunay
gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)
gmsh.option.setNumber('Mesh.Format', 1)
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 50)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.model.mesh.generate(3)
#gmsh.fltk.run()
# mesh = 3d topology
# ct = volume tags
# fm = surface tags
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, modelRank, gdim=3)
if modelRank == mpiRank:
gmsh.finalize()
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
elem = ufl.FiniteElement('Nedelec 1st kind H(curl)', mesh.ufl_cell(), degree=2)
V = fem.FunctionSpace(mesh, elem)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Print out BC tags
with io.XDMFFile(mesh.comm, "BCs.xdmf", "w") as xx:
xx.write_mesh(mesh)
xx.write_meshtags(fm)
n = ufl.FacetNormal(mesh)
# Dirichlet BCs
facets = fm.find(3) # PEC facets
ubc = fem.Function(V)
ubc.x.set(0+0j)
dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
bc = fem.dirichletbc(ubc, dofs)
# Set up forms for eigenvalue problem Ax -k0^2*Bx = 0
a = fem.form(ufl.inner(ufl.curl(u), ufl.curl(v)) * ufl.dx)
b = fem.form(ufl.inner(u, v) * ufl.dx)
A = fem.petsc.assemble_matrix(a, bcs=[bc])
A.assemble()
B = fem.petsc.assemble_matrix(b, bcs=[bc])
B.assemble()
eigs = slepc.EPS().create(mesh.comm)
eigs.setOperators(A, B)
eigs.setProblemType(slepc.EPS.ProblemType.GHEP) # General hermitian
eigs.setTolerances(1.0e-8)
eigs.setType(slepc.EPS.Type.KRYLOVSCHUR)
st = eigs.getST() # Invoke spectral transformation
st.setType(slepc.ST.Type.SINVERT)
eigs.setWhichEigenpairs(slepc.EPS.Which.TARGET_REAL) # target real part of eigenvalue
eigs.setTarget(-1.0)
eigs.setDimensions(nev=5) # Number of eigenvalues
eigs.solve()
eigs.view()
eigs.errorView()
sys.exit(0)
