Another question about FenicsX boundary conditions

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?

A general question:

  • Why aren’t you using gmsh to thg the boundaries, as you are using gmsh to generate the mesh?

Then to the issue:

  • When you use locate_entities_boundary, it checks all vertices of that the logical function you sent in evaluates to true or false. If all are true, it returns the entity. If not, it does not return it. What I would do if I were you would be the following:
  1. Create a list of all facets in the mesh, i.e. facets = np.arange(mesh.topology.index_map(mesh.topology.dim-1).size_local + mesh.topology.index_map(mesh.topology.dim-1).num_ghosts, dtype=np.int32) and a list of values `values=np.zeros_like(facets, dtype=np.int32)
  2. Create a list of all facets on the boundary, i.e. use bndy_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology), and then call values[bndy_facets] = 1
  3. Use locate_entities_boundary for the tag 2 and 3 and modify the values as done in step 2.

This is only if you really do not want to use gmsh for this tagging.

1 Like

Thanks! I took your suggestion to apply the boundary tags in GMSH and it works great.

The next problem is getting the boundary integrals to work properly. I defined a vector function Hinc in a class and attempt to perform a boundary integral over the appropriate surface (with boundary tag β€œ2”). Code follows:

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
    k0 = 1.92
    eta0 = 377.0
    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)
# Set up BC tags
# Boundary conditions: 1=PEC, 2=input and 3=output BCs
       Shield, InPort, OutPort = 1, 2, 3
       bb = gmsh.model.getEntities(dim=2)
       print(bb)
       pec = []
       inportA = []
       outportA = []
       bndry = gmsh.model.getBoundary(bb, oriented=False)
       for bnd in bb:
          CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
          if np.allclose(CoM, [a/2, -lc, w]):
             inportA.append(bnd[1])
          elif np.allclose(CoM, [a/2, b/2, l]):
             outportA.append(bnd[1])
          else:
             pec.append(bnd[1])
       print(pec)
       print(inportA)
       print(outportA)
       gmsh.model.addPhysicalGroup(2, pec, Shield)
       gmsh.model.setPhysicalName(2, Shield, "PEC")
       gmsh.model.addPhysicalGroup(2, inportA, InPort)
       gmsh.model.setPhysicalName(2, InPort, "Inlet")
       gmsh.model.addPhysicalGroup(2, outportA, OutPort)
       gmsh.model.setPhysicalName(2, OutPort, "Outlet")
       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)

    
    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(fm)
    print(PETSc.ScalarType)
#Set up weak form
    L = (ufl.inner(ufl.curl(u), ufl.curl(v)) - k0 * k0 * Dk * ufl.inner(u, v)) * ufl.dx

# Facet normal unit vectors
    n = ufl.FacetNormal(mesh)
# Surface subdomains
    ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
#Import BCs into system
# Incident field
    class Hinc:
        def __init__(self, a: float, w: float, H: complex):
            self.a = a
            self.w = w
            self.H = H
        def eval( \
                self, x: np.typing.NDArray[np.float64]) -> tuple[ \
                    np.typing.NDArray[np.complex128],\
                    np.typing.NDArray[np.complex128],\
                    np.typing.NDArray[np.complex128]]:
                hx = H * (x[2]-w)/(2.0*np.pi*np.sqrt(np.pow(x[0]-0.5*a, 2.0) + np.pow(x[2]-w, 2.0)))
                hy = 0.0
                hz = H * (x[0]-0.5*a) / (2.0 * np.pi*np.sqrt(np.pow(x[0]-0.5*a,2.0)+np.pow(x[2]-w,2.0)))
                return(hx, hy, hz)

#Build RHS
    uh = fem.Function(V, dtype=np.complex128)
    ac = fem.Constant(mesh, PETSc.ScalarType(a))
    wc = fem.Constant(mesh, PETSc.ScalarType(w))
    f = Hinc(ac, wc, 1j)
    uh.interpolate(f.eval)
    rhs = 2.0 * k0 * eta0 * inner(cross(n, uh), v) * ds(2)

# Dirichlet n x E = 0
    facets = fm.find(1)
    dofs = fem.locate_dofs_topological(V, mesh.topology.dim-1, facets)
    bc = fem.dirichletbc(np.array((0,0,0), dtype=PETSc.ScalarType), dofs)

    return (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)

I get the following error:

Antenna height= 0.648000, Antenna distance from short= 0.554000
[(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9), (0, 10), (0, 11), (0, 12), (0, 13)]
[(0, 11), (1, 16), (2, 9)]
[(0, 7), (0, 8), (1, 8), (1, 9), (2, 3)]
[(0, 12), (0, 13), (1, 18), (1, 19), (2, 11)]
[(2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7), (2, 8), (2, 9), (2, 10), (2, 11), (2, 12)]
[1, 2, 3, 5, 6, 7, 8, 9, 10, 12]
[11]
[4]
<class 'numpy.complex128'>
Traceback (most recent call last):
  File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 343, in interpolate
    _interpolate(u, cells)
  File "/usr/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 319, in _interpolate
    self._cpp_object.interpolate(u, cells)
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfinx.cpp.fem.Function_complex128, f: numpy.ndarray[numpy.complex128], cells: numpy.ndarray[numpy.int32]) -> None
    2. (self: dolfinx.cpp.fem.Function_complex128, u: dolfinx.cpp.fem.Function_complex128, cells: numpy.ndarray[numpy.int32]) -> None
    3. (self: dolfinx.cpp.fem.Function_complex128, expr: dolfinx::fem::Expression<std::complex<double> >, cells: numpy.ndarray[numpy.int32]) -> None

Invoked with: <dolfinx.cpp.fem.Function_complex128 object at 0x7ff470612e30>, <bound method Cost.<locals>.Hinc.eval of <__main__.Cost.<locals>.Hinc object at 0x7ff46fc0b1c0>>, array([    0,     1,     2, ..., 37424, 37425, 37426], dtype=int32)

Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
when compiling your pybind11 module.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/bill/Cluster/Fenics2020/CoaxWGCoupler/CoaxWGTrans.py", line 183, in <module>
    res = Cost(x0)
  File "/home/bill/Cluster/Fenics2020/CoaxWGCoupler/CoaxWGTrans.py", line 162, in Cost
    uh.interpolate(f.eval)
  File "/usr/lib/petscdir/petsc3.15/x86_64-linux-gnu-complex/lib/python3/dist-packages/dolfinx/fem/function.py", line 348, in interpolate
    self.interpolate(u(x), cells)
  File "/home/bill/Cluster/Fenics2020/CoaxWGCoupler/CoaxWGTrans.py", line 152, in eval
    hx = H * (x[2]-w)/(2.0*np.pi*np.sqrt(np.pow(x[0]-0.5*a, 2.0) + np.pow(x[2]-w, 2.0)))
NameError: name 'H' is not defined
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[13100,1],0]
  Exit code:    1
--------------------------------------------------------------------------

I am either not defining the function correctly or am not doing the eval correctly in the boundary integration. H is indeed defined as β€œ1j”. Something is probably fishy with the variable passing. What is the best way to do this?

H is not defined.
self.H however, is defined. This is how classes in python work, see for instance: self in Python class - GeeksforGeeks

Also note that numpy does not implement pow as you are calling it.
instead of calling

you should call (x[0]-0.5*a)**2 and similar calls.

Finally, you are not returning the correct type for hy. It should be a numpy array, for instance

should be

            return (hx, np.full_like(hx, hy, dtype=np.complex128), hz)

There are also issues in your DirichletBC assignment (it does not relate to a function space:

and should be

    ubc = fem.Function(V)
    ubc.x.set(0+0j)
    bc = fem.dirichletbc(ubc, dofs)

1 Like

Thanks for that. Sometimes, I miss some of the subtleties of the layers of abstraction in Python. Anyway, the code below now runs without errors, but the solution is wrong. There is still an error in the invocation of the β€œHinc” boundary condition. The E field should be azimuthal about the coaxial input port and it is not. I think I will take another look at it tomorrow.

EDIT 2022-01-17: Got it all to work (working code below). Thanks to Dokken for his assistance.

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 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]
    lcc = 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
    eps_0 = 1.0
    eps_c = 2.2
    k0 = 2.09
    beta = np.sqrt(k0 * k0 - (np.pi / a)**2)
    eta0 = 377.0
    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, lcc)
# Set up BC tags
# Boundary conditions: 1=PEC, 2=input and 3=output BCs
       Shield, InPort, OutPort = 1, 2, 3
       bb = gmsh.model.getEntities(dim=2)
       print(bb)
       pec = []
       inportA = []
       outportA = []
       bndry = gmsh.model.getBoundary(bb, combined=True, oriented=False)
       print(bb)
       print(bndry)
       for bnd in bb:
          CoM = gmsh.model.occ.getCenterOfMass(bnd[0], bnd[1])
          print(CoM)
          if np.allclose(CoM, [a/2, -lc, w]):
             inportA.append(bnd[1])
          elif np.allclose(CoM, [a/2, b/2, l]):
             outportA.append(bnd[1])
          elif not np.allclose(CoM, [a/2, 0, w]):
             pec.append(bnd[1])
       print(pec)
       print(inportA)
       print(outportA)
       gmsh.model.addPhysicalGroup(2, pec, Shield)
       gmsh.model.setPhysicalName(2, Shield, "PEC")
       gmsh.model.addPhysicalGroup(2, inportA, InPort)
       gmsh.model.setPhysicalName(2, InPort, "Inlet")
       gmsh.model.addPhysicalGroup(2, outportA, OutPort)
       gmsh.model.setPhysicalName(2, OutPort, "Outlet")
       gmsh.model.occ.synchronize()
       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)
    print(PETSc.ScalarType)

    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, eps_c, dtype=PETSc.ScalarType)
    Dk.x.array[Air] = np.full_like(Air, eps_0, dtype=PETSc.ScalarType)

    
    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(fm)

# Facet normal unit vectors
    n = ufl.FacetNormal(mesh)
# Surface subdomains
    ds = ufl.Measure("ds", domain=mesh, subdomain_data=fm)
#Import BCs into system
# Incident field
    class Hinc:
        def __init__(self, a, w, H):
            self.a = a
            self.w = w
            self.H = H
        def eval(self, x):
            hx = -self.H * (x[2]-self.w)/(2.0*np.pi*np.sqrt((x[0]-0.5*self.a)**2.0 + (x[2]-self.w)**2.0))
            hy = np.full_like(hx, 0.0+0j, dtype=np.complex128)
            hz = self.H * (x[0]-0.5*self.a) / (2.0 * np.pi*np.sqrt((x[0]-0.5*self.a)**2.0+(x[2]-self.w)**2.0))
            return(hx, hy, hz)

# Dirichlet n x E = 0
    facets = fm.find(1)
    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)

#Build RHS
    uh = fem.Function(V, dtype=np.complex128)
    f = Hinc(a, w, 1+0j)
    uh.interpolate(f.eval)
    f_inc = 2.0 * k0 * eta0 * ufl.inner(ufl.cross(n, uh), v) * ds(2)
    uh.name = "Hinc"
    with io.XDMFFile(mesh.comm, "Hinc.xdmf", "w") as xdmf:
       xdmf.write_mesh(mesh)
       xdmf.write_function(uh)

#Set up weak form
    L = (ufl.inner(ufl.curl(u), ufl.curl(v)) - k0 * k0 * Dk * ufl.inner(u, v)) * ufl.dx 

# Build outgoing wave boundaries
# Input port
    bc1 = 1j * k0 * np.sqrt(eps_c) * ufl.inner(ufl.cross(n, u), ufl.cross(n,v)) * ds(2)
    L += bc1
# Output port
    bc2 = 1j * beta * ufl.inner(ufl.cross(n, u), ufl.cross(n,v)) * ds(3)
    L += bc2
# Solve system
    Lin_system = fem.petsc.LinearProblem(L, f_inc, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    E = Lin_system.solve()
    lu_solver = Lin_system.solver
    lu_solver.view()

    E.name = "ElectricField"
    with io.XDMFFile(mesh.comm, "Er.xdmf", "w") as xdmf:
       xdmf.write_mesh(mesh)
       xdmf.write_function(E)

# Print out time steps
    Et = fem.Function(V)
    Et.name = "AnimEfield"

    with io.XDMFFile(mesh.comm, "Movie.xdmf", "w") as xx:
        xx.write_mesh(mesh)
        for t1 in range(50):
           Et.vector.array = E.vector.array * np.exp(1j * np.pi * t1 /25.0)
           xx.write_function(Et, t1)

    return (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)