Bug with 3D prism edges

[DerekK44]

Hello,

I'd like to report what seems to be a bug I've noticed for the subpixel averaging of outer edges of prisms with large inner surfaces. For example, a prism with a narrow opening and large inner surface (such as the shape drawn below) will have the issue:

The issue seems to be that when meep (or libctl) builds the structure or does the subpixel averaging, some regions of the polygon have a shorter dimension than others, as shown by looking at slices:

This issue becomes especially problematic when the prism thickness is on the order of a few pixels. I have attached a MWE to help reproduce the issue above, as well as shows another example with issues generating the following polygon (of width 2 pixels large):

Let me know if I can provide any other information. Sorry I can't be of more help in finding the source of the issue. I'm still looking through meep/libctl to understand how things work.

[stevengj]

cc @HomerReid

[m9631m]

Just a suggestion, sometimes Paraview interpolates the points to make it nicer. What you are saying, is clearly an issue but for 2D outputs, I usually plot it with points (with large point size) to make sure it is not visualization issue.

[DerekK44]

Thanks for the suggestion, m9631m. Attached are the 3D images using the "points" rather than a "surface". We can see the features I reported above in both visualization modes:

[HomerReid]

Thanks for the well-documented bug report and especially for being the first to do MEEP calculations involving recognizers from TRON, the greatest sci-fi film of all time.

Two questions:
(1) How did you make those nice plots in paraview, using points and using surfaces?
(2) Your prisms are extruded in the Y direction. Does the issue remain if you rotate so that the prisms are extruded in the Z direction?

For the time being, I haven't exactly reproduced your result, but I do observe a discrepancy in the epsilon grids with and without averaging. The following images show constant-Y slices at two heights: one at the lower boundary of the object, and one slicing through the object.

Lower boundary:

Slicing through object:

[stevengj]

Thanks Homer, that's informative. It looks like there is a bug in the point_in_object routine for the prism.

[HomerReid]

In libctl#18 I did find and fix a subtle issue in the point-in-prism algorithm that was responsible for the spurious narrow strips of purple in the figures above.

However, I'm not sure whether or not this fixes the original issue. Here's a threshold image, similar to the above, produced after fixing the above issue.

It looks like the issue may still be there.

When you get a chance, could you double-check whether or not the updated version of libctl makes any difference on your geometries? It would also be useful to know if anything changes when the prism axis is the z-axis rather than the y-axis.

[DerekK44]

Hi Homer,

Thanks for looking into this. To answer your questions above:

1. First, I converted the epsilon .h5 to .vtk (using h5tovtk), then opened the .vtk in paraview. Then, with the .vtk file selected in the "Pipeline Browser" window, I applied a "Threshold". Under the Threshold options, I set a range (max and min of the data to be shown). This range is what I included in the text annotations above.

2. I checked both y-extruded and z-extruded prisms before and after rebuilding with the changes you made to libctl. Unfortunately the problem seems to persist in all the cases I tested.

[DerekK44]

I double checked the math in intersect_line_with_segment, and it seems right to me..

Is it possible there's a subtle bug somewhere in the intersect_line_with_prism or normal_to_prism routines of geom.c? Especially since the problem occurs at the surfaces of the prism, including the top/bottom prism surfaces. I'll keep looking

[HomerReid]

Yes, it seems likely that your geometry exposes a bug in one of the new prism geometry routines in libctl (specifically in utils/geom.c) that wasn't hit by previous geometries. My strategy for pursuing it has been to instantiate your structure in two ways, as a prism and as a union of blocks, and look for differences in the results of calls to normal_to_object and other routines. Let me know if you come across any further insights.

[DerekK44]

Ok, I think I have an idea of the problem.

In normal_to_prism (lines 2339 - 2369) the normal distance between each prism surface and a point is computed. It seems to assume the plane of the prism surface extends to infinity, which is where I think the problem is.

As a result, a prism surface that is far away from the point can be chosen as the "closest surface" if the normal of the infinite plane is sufficiently small. See image attached. For point xp that lies on the top surface (in the out-of-page direction) of a 3D 'recognizer' :) the top surface of the prism should be chosen, and the normal-direction should be out-of-page. However, the surface in blue could actually be selected (even though it is far away, the normal-distance from the plane is zero).

Perhaps the normal distance should be used only if the normal vector (containing the point) intersects with the plane on the true surface (i.e. the blue region shown above). However, if the normal vector does not intersect with the true surface, then the shortest Euclidean distance from the edges (like points v1 and v2 above) of the true surface could be used?

[HomerReid]

Many thanks for that brilliant piece of detective work! It would have taken me a very long time to identify this issue. I implemented a fix in https://github.com/stevengj/libctl/pull/19 that does seem to correct the problem.

[DerekK44]

Thanks for the fix! I rebuilt libctl using your specific branch and do get the same result for the MWE using -case=0. Could you see if -case=1 and -case=2 are also fixed? When I tried I seem to still be getting the same bug for those cases..

[HomerReid]

I've updated libctl#19 with a revised implementation of normal_to_object for prisms. It now seems to be doing the right thing on your case-2 geometry:

but I'm not sure if it's doing the right thing on the recognizers.

[HomerReid]

Just FYI, the updates to point_in_prism and normal_to_object for prisms have now been merged into libctl. Here are pictures of the recognizer with resolution={30,60}:

• Resolution 30:

• Resolution 60: