[BUG]: (Fixed in this post) The "eps" function fails when the "Omega" functions are tweaked in the tutorial "Electromagnetic modal analysis for a waveguide"

[FyDininno]

Summarize the issue

Changing the domain conditions of the "Omega" functions just a little bit completely breaks the simulation. Instead of using "Omega" functions to assign the values to eps, what you have to do is write an "eps_expr(x)" function which returns the scalar value of eps over the domain at each x. Next, you have to run "eps.interpolate(eps_expr)" to correctly assign the dielectric values to the eps function. This way, it does not break.

How to reproduce the bug

To reproduce the bug, use these omega functions instead of the ones given in the tutorial, and set all of the dielectric permittivities to 1 just to illustrate that the bug arises from the assignment of eps, and not the values of the permittivities. As you can see, simply by making the dielectric domain two dimensional instead of one dimensional completely breaks the simulation, giving nonsensical results.

Minimal Example (Python)

eps_v = 1
eps_d = 1

def Omega_d(x): # Returns True within the correct domain, false outside the domain
    return (x[0] <= w/2) & (x[1] <= h/2)

def Omega_v(x):
    return (x[0] >= w/2) | (x[1] >= h/2)

D = fem.functionspace(msh, ("DQ", 0)) # This defines the FEM function space over the mesh with linear (degree 0) discontinuous functions defined on quadilateral elements
eps = fem.Function(D) # A template function created by elements within D, as of now, all the basis functions have a coefficient of zero

cells_v = locate_entities(msh, msh.topology.dim, Omega_v)
cells_d = locate_entities(msh, msh.topology.dim, Omega_d) # locates the cells which have the correct truth value

eps.x.array[cells_v] = np.full_like(cells_v, eps_v, dtype=scalar_type) # .x gives acces to the DOF vector for eps. The argument specifies the cell values that are being modified. The numpy array formats an array to the likes of the argument with scalar values of the correct dielectric permittivities in their place.
eps.x.array[cells_d] = np.full_like(cells_d, eps_d, dtype=scalar_type)

Output (Python)

No response

Version

main branch

DOLFINx git commit

No response

Installation

I used conda and set the vscode environment to use the fenicsx conda environment.

Additional information

If you simply add

def eps_expr(x):
  cond_d = Omega_d(x)
return np.where(cond_d, eps_d, eps_v)

D = fem.functionspace(msh, ("CG", 1))
eps = fem.Function(D)
eps.interpolate(eps_expr)

It should work