I have finally decided to take the plunge and migrate from old FEnICS to new FEnICSx. One area that is very different is the invocation of boundary conditions. I am almost there, bit there is a problem that has stumped me.
I am migrating one of my old models (Variation of coax-to-waveguide model found here) to FEnICSx and the boundary facet labeling scheme is not behaving as I expect.
First, the new FEnICSx code:
import sys
import numpy as np
from scipy import optimize as opt
from mpi4py import MPI as nMPI
from petsc4py import PETSc
import gmsh
import meshio
import ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx import fem
def Cost(x):
#Set up processor rank for job control
comm = nMPI.COMM_WORLD
mpiRank = comm.rank
model_rank = 0
h = x[0]
w = x[1]
lc = 0.05
lm = 0.2
ls = 0.02
a = 2.286 # WG width
b = 1.016 # WG height
l = 4.0 # WG length
r = 0.065 # probe radius
rc = 0.22 # Coax outer rad
lc = 1.0 # Coax length
eps = 1.0e-3
if mpiRank == model_rank:
print("Antenna height= {0:<f}, Antenna distance from short= {1:<f}".format(h, w))
gmsh.initialize()
gmsh.option.setNumber('General.Terminal', 0)
gmsh.model.add("WGLaunch")
gmsh.model.setCurrent("WGLaunch")
gmsh.model.occ.addBox(0.0, 0.0, 0.0, a, b, l, 1) # WG shield
gmsh.model.occ.addCylinder(a/2, -lc, w, 0.0, lc, 0.0, rc, 2) # Coax shield
gmsh.model.occ.addCylinder(a/2, -lc, w, 0.0, lc+h, 0.0, r, 3) # Coax center cond & probe
gmsh.model.occ.cut([(3,1)],[(3,3)], 4, removeObject=True, removeTool=False)
gmsh.model.occ.cut([(3,2)],[(3,3)], 5, removeObject=True, removeTool=True)
gmsh.option.setNumber('Mesh.MeshSizeMin', ls)
gmsh.option.setNumber('Mesh.MeshSizeMax', lm)
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.Smoothing', 100)
gmsh.option.setNumber('Mesh.MinimumCirclePoints', 36)
gmsh.option.setNumber('Mesh.CharacteristicLengthFromCurvature', 1)
gmsh.model.occ.fragment([(3,4),(3,5)],[], -1)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(3, [4], 1)
gmsh.model.addPhysicalGroup(3, [5], 2)
gmsh.model.setPhysicalName(3, 1, "Air")
gmsh.model.setPhysicalName(3, 2, "Diel")
pt = gmsh.model.getEntities(0)
print(pt)
gmsh.model.mesh.setSize(pt, lm)
ov = gmsh.model.getEntitiesInBoundingBox(a/2-2*r-eps, h-eps, w-2*r-eps, a/2+2*r+eps, h+eps, w+2*r+eps)
print(ov)
gmsh.model.mesh.setSize(ov, ls)
sv = gmsh.model.getEntitiesInBoundingBox(a/2-rc-eps, -eps, w-rc-eps, a/2+rc+eps, eps, w+rc+eps)
print(sv)
gmsh.model.mesh.setSize(sv, ls)
bv = gmsh.model.getEntitiesInBoundingBox(a/2-rc-eps, -lc-eps, w-rc-eps, a/2+rc+eps, -lc+eps, w+rc+eps)
print(bv)
gmsh.model.mesh.setSize(bv, lc)
gmsh.model.mesh.generate(3)
gmsh.model.mesh.optimize("Gmsh")
# gmsh.write("CoaxWG.msh")
# gmsh.fltk.run()
mesh, ct, fm = gmshio.model_to_mesh(gmsh.model, comm, model_rank, gdim=3)
if mpiRank == model_rank:
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)
# Material type
Q = fem.FunctionSpace(mesh, ("DG", 0))
Dk = fem.Function(Q)
Diel = ct.find(2)
Air = ct.find(1)
Dk.x.array[Diel] = np.full_like(Diel, 2.2, dtype=PETSc.ScalarType)
Dk.x.array[Air] = np.full_like(Air, 1.0, dtype=PETSc.ScalarType)
# Boundary conditions: 1=PEC, 2=input and 3=output BCs
boundaries = [(1, lambda x: np.logical_not(np.logical_or(np.isclose(x[1], -lc, rtol=1.0e-6), np.isclose(x[2], l, rtol=1.0e-6)))), \
(2, lambda x: np.isclose(x[1], -lc)), \
(3, lambda x: np.isclose(x[2], l))]
facet_indices = []
facet_markers = []
for (marker, locator) in boundaries:
facets = locate_entities_boundary(mesh, mesh.topology.dim-1, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, mesh.topology.dim-1, facet_indices[sorted_facets], facet_markers[sorted_facets])
from dolfinx import io
with io.XDMFFile(mesh.comm, "Dk.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(Dk)
with io.XDMFFile(mesh.comm, "BCs.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_meshtags(facet_tag)
return (0.0)
#h1 = 0.648+0.14 # Starting probe height
#w1 = 0.554-0.25 # Starting probe distance from short
h1 = 0.648 # Starting probe height
w1 = 0.554 # Starting probe distance from short
#Set up optimization
#bounds = opt.Bounds([h1-0.3, w1-0.4], [h1+0.3, w1+0.4]) # The trial range
x0 = np.array([h1, w1]) # Starting position deviation
#res = opt.minimize(Cost, x0, method='trust-constr', options={'verbose': 1}, bounds=bounds)
#print(res)
res = Cost(x0)
sys.exit(0)
What I want to do is label an input boundary (at x[1] == -1) and an output boundary (at x[2] == 4) and all the rest and a Dirichlet condition that is set to zero. The preceding code finds the input and output boundaries (with markers 2 and 3). However, the “all the rest” (facets not assigned to marker 2 or 3) should be labelled “1”, but from the Paraview threshold plot, is seems that not all the facets that should lie in “1” (the blue surface) near the “2” and “3” surfaces are labeled.
Facets with marker “1” do not seem to be all there. There is a gap between two surfaces. Is this a plotting artifact or have the facets that are not seen actually been left out? How do I fix this?
