Prism geometry rendering artifacts

[ianwilliamson]

I am encountering some strange artifacts in Meep's rendering of Prism geometries onto the simulation grid. Consider the case shown below. The prism vertices are traced on the left and the Meep simulation domain is shown on the right, with the same vertices overlaid as x's.

You can see that there is an artifact around the lower left part of the cutout. This seems to be coming from the extra vertex that sits on that diagonal interface. While this vertex is not strictly necessary to describe the outline of this prism, I would not have expected Meep to render this geometry as shown in the image.

The code for reproducing the result shown above is provided below:

import meep as mp
import matplotlib.pyplot as plt
geometry = [
    mp.Prism(
        vertices=[
            mp.Vector3(3450.0008325358162, 300.0, 0.0),
            mp.Vector3(3450.0008209110665, 235.0, 0.0),
            mp.Vector3(3450.000676495217, 195.0, 0.0),
            mp.Vector3(3240.0000608065666, 115.0, 0.0),
            mp.Vector3(3230.0000572294366, 54.999999999999545, 0.0),
            mp.Vector3(3245.0, 39.99988881841773, 0.0),
            mp.Vector3(3255.0, 30.000042326215862, 0.0),
            mp.Vector3(3450.00077515856, 125.0, 0.0),
            mp.Vector3(3450.0004644694577, -100.0, 0.0),
            mp.Vector3(3200.0, -100.0, 0.0),
            mp.Vector3(3200.0, 300.0, 0.0),
        ],
        height=100,
        sidewall_angle=0,
        material=mp.Medium(index=4),
    ),
]

simulation = mp.Simulation(
    cell_size=mp.Vector3(300.00099904297946, 480.0, 0.0),
    boundary_layers=[],
    geometry=geometry,
    geometry_center=mp.Vector3(3322.401246760824, 95.34725790541609, 0.0),
    resolution=0.1,
    eps_averaging=False,
)
simulation.init_sim()
epsr = simulation.get_array(
    center=simulation.geometry_center,
    size=simulation.cell_size,
    component=mp.Dielectric,
    cmplx=False,
)

(x, y, z, _) = simulation.get_array_metadata(
    center=simulation.geometry_center,
    size=simulation.cell_size,
)

fig, ax = plt.subplots(1,2, figsize=(6,6), constrained_layout=True)
ax[1].pcolormesh(x, y, epsr.T)

for g in simulation.geometry:
  if isinstance(g, mp.Prism):
    ax[1].plot([vtx.x for vtx in g.vertices], [vtx.y for vtx in g.vertices], 'rx')
    ax[0].plot([vtx.x for vtx in g.vertices], [vtx.y for vtx in g.vertices], 'r-')

ax[0].axis('image')
ax[0].set_title('Polygon')
ax[1].axis('image')
ax[1].set_title('Meep simulation domain')

[oskooi]

The artifact seems to be caused by floating-point rounding error in the point_in_prism function of libctl since applying a small perturbation to the vertex at mp.Vector3(3245.0, 39.99988881841773, 0.0) (the "unnecessary" vertex on the diagonal) to change it to e.g. mp.Vector3(3245.0, 40.0, 0.0) produces the correct geometry without any artifacts.

[ianwilliamson]

Thanks for looking into this! Applying the following function:

import numpy as onp

def round(prism, decimals=4):
  rounded_vtxs = []
  for vtx in prism.vertices:
    rounded_vtxs.append(mp.Vector3(*onp.around(vtx, decimals)))
  return mp.Prism(
      vertices=rounded_vtxs,
      height=prism.height,
      material=prism.material,
  )

to the geometry:

geometry = [round(g, decimals=2) for g in geometry]

does seem to eliminate the artifact:

Perhaps the calculation in point_in_prism has a bug if it is unstable to floating point rounding errors?

[smartalecH]

The point in polygon problem (which is how point_in_prism works) is a notoriously hard one to solve (i.e. there's no one-size-fits-all solution). See #864 for some background.

Long story short, if your polygon is too complicated (lots of vertices, several concave sides, etc.) it may fail to properly render.

There are some very robust (and slow) libraries that can handle these corner cases well, but often it's easier to do a little geometry preprocessing.

[ianwilliamson]

The point in polygon problem (which is how point_in_prism works) is a notoriously hard one to solve (i.e. there's no one-size-fits-all solution). See #864 for some background.

Is the example that I've shown above an edge case though? The last message of the issue you linked seems to indicate that a floating point round off bug still persists. That seems to be what's going on here.

[smartalecH]

Sorry, I should clarify. In addition to my previous comment, I don't think our current algorithm is backwards stable which is why it's so sensitive to round off error and floating point precision. By edge cases, I mean anything that triggers this sensitivity (similar to Gaussian Elimination or Gram-Schmidt).

As I mentioned, this is a common issue with this particular problem. Mainstream libraries, like CGAL, get around this by doing lots of checking. Of course, this induces a lot of slow branching.

In other words, I don't think this is a bug, but rather a consequence of the algorithm.

But I'd be happy if I were wrong!

[oskooi]

Fixed by NanoComp/libctl#59?

[ahoenselaar]

No. That just prevents segfaults under certain circumstances.

…

On Tue, Nov 9, 2021 at 8:12 PM Ardavan Oskooi ***@***.***> wrote: Fixed by NanoComp/libctl#59 <NanoComp/libctl#59>? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1523 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA4VK535AH26ZQ2L2Z6J4JLULHWJHANCNFSM4Y4LO74A> .