Hi all,
When setting up an eigenvalue problem that uses periodic boundary conditions, I am having trouble assembling the matrix eigenvalue problem using MPC. On matrix assembly:
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 dolfinx_mpc
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
BCweight = 100.0
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)
V = fem.VectorFunctionSpace(mesh, ("CG", 2))
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)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
# Periodic BC mapping
def PeriodicMapping(x):
out = np.zeros(x.shape)
out[0] = x[0]
out[1] = x[1]
out[2] = x[2] + l
return out
# Dirichlet BCs using weak invokation
Lb = BCweight * ufl.inner(ufl.cross(n, u), ufl.cross(n, v)) * ds(3)
bc=[]
# Create multipoint boundary. Master boundary has tag 2,
# slave has tag 1 - set on periodic constraint.
# (see GMSH physical group definitions above)
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_topological(V, fm, 1, PeriodicMapping, bc, 1.0)
mpc.finalize()
# Set up forms for eigenvalue problem Ax -k0^2*Bx = 0
az = ufl.as_vector((0, 0, 1)) # Z-directed unit vector
a = fem.form(ufl.inner(Gamma * ufl.cross(az, u) + ufl.curl(u), Gamma * ufl.cross(az, v) + ufl.curl(v)) * ufl.dx + Lb)
b = fem.form(ufl.inner(u, v) * ufl.dx)
num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
num_dofs_global = V.dofmap.index_map.size_global * V.dofmap.index_map_bs
print(f"Number of dofs (owned) by rank {mpiRank}: {num_dofs_local}")
if mpiRank == modelRank:
print(f"Number of dofs global: {num_dofs_global}")
A = dolfinx_mpc.assemble_matrix(a, mpc, bcs=bc, diagval=1.0)
B = dolfinx_mpc.assemble_matrix(b, mpc, bcs=bc, diagval=1.0)
eigs = slepc.EPS().create(mesh.comm)
eigs.setOperators(A, B)
eigs.setProblemType(slepc.EPS.ProblemType.GNHEP) # General non-hermitian
eigs.setTolerances(1.0e-8)
eigs.setType(slepc.EPS.Type.KRYLOVSCHUR)
st = eigs.getST() # Invoke spetcral 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)
I get the following error:
File "/home/bill/Cluster/Fenics2020/Helix/HelicalSection3.py", line 143, in <module>
A = dolfinx_mpc.assemble_matrix(a, mpc, bcs=bc, diagval=1.0)
File "/usr/lib/python3/dist-packages/dolfinx_mpc/assemble_matrix.py", line 51, in assemble_matrix
A = cpp.mpc.create_matrix(form, constraint[0]._cpp_object,
TypeError: create_matrix(): incompatible function arguments. The following argument types are supported:
1. (arg0: dolfinx.cpp.fem.Form_float64, arg1: dolfinx_mpc.cpp.mpc.MultiPointConstraint) -> mat
2. (arg0: dolfinx.cpp.fem.Form_float64, arg1: dolfinx_mpc.cpp.mpc.MultiPointConstraint, arg2: dolfinx_mpc.cpp.mpc.MultiPointConstraint) -> mat
It seems I am doing something wrong in invoking the MPC. What am I missing?