Hello all,
I am presently attempting to solve the following weak form in Dolfinx (which I have successfully solved using legacy Dolfin, which is discussed here):
F=\left (\Gamma\hat{\mathbf{z}}\times\mathbf{E} + \nabla\times\mathbf{E}\right )^*\cdot\left (\Gamma\hat{\mathbf{z}}\times\mathbf{E}+\nabla\times\mathbf{E}\right ) - k_0^2\mathbf{E}^*\cdot\mathbf{E},
where the periodic boundary is enforced between z=0 and z=l. The term \Gamma is the complex valued Floquet propagation factor (for lossless structures, we assume to take imaginary values only).
I have a couple of problems -
- In line 132 of the code (where I define the MPC constraint), there is a problem with the variable types, but after looking at the man page for “create_periodic_constraint_topological” and looking at examples, something is still not right. Perhaps my boundary relation function output is not correctly presented.
- What is the best way to set up the resulting eigenvalue problem?
# Generate unit cell of helical structure
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
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)
# 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
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)
# 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.asvector((0, 0, 1)) # Z-directed unit vector
a = ufl.inner(Gamma * ufl.cross(az, u) + ufl.curl(u), Gamma * ufl.cross(az, v) + ufl.curl(v)) * dx
b = ufl.inner(u, v) * dx
A = dolfinx_mpc.assemble_matrix(a, mpc, bc, 1.0)
A.assemble()
B = dolfinx_mpc.assemble_matrix(b, mpc, bc, 1.0)
B.assemble()
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)