Bug with 3D prism edges #400
Comments
|
cc @HomerReid |
|
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. |
|
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:
|
|
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: 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:
|
|
Thanks Homer, that's informative. It looks like there is a bug in the |
|
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. |
|
Hi Homer, Thanks for looking into this. To answer your questions above:
|
|
I double checked the math in Is it possible there's a subtle bug somewhere in the |
|
Yes, it seems likely that your geometry exposes a bug in one of the new prism geometry routines in libctl (specifically in |
|
Ok, I think I have an idea of the problem. In 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
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? |
DerekK44
changed the title
|
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.
|
|
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.. |
|
I've updated libctl#19 with a revised implementation of
but I'm not sure if it's doing the right thing on the recognizers. |
|
Just FYI, the updates to
|
|
Hi Homer, Thanks for the updates. I think in the -case=2 image above, there's an area actually missing (below they are side-by-side): I found one potential issue in the
One potential fix is to first check if the line goes through a vertex (i.e. if I am not sure if this is the bug that is causing all the weird effects on the recognizers, but this seemed to cause issues when I rotate the recognizer by 45 degrees. I'll keep looking and see if I find anything else! |
|
Ok, so I suspect the second bug (that is causing the edges to be off in the 'recognizers') is happening in the Might this be the cause of the above issue? |
|
Thanks again for all the detective work and the extremely useful bug reports. The issue with the plumb lines intersecting corners is subtle indeed. Thanks for identifying it, and for proposing a fix. I will patch the code to take this into account, but I suspect it's not the primary cause of the main issue you reported---it seems like it would only affect a small subset of points that happen to satisfy this geometrical coincidence, whereas the issue you report seems to extend over wider regions. The second possibility you raise is also a plausible culprit. The idea of retaining only the nearest two intersections was proposed by @stevengj (see e.g. here), with the idea being that subpixel averaging is entirely local and cannot be affected by distant surface intersections. However, it seems possible that this reasoning only applies for star-shaped objects, or at any rate is invalid for objects with inner surfaces like in this case. One observation is that the issue seems to be present for the recognizers only when the claw feet are present, and not for the simple inverted-U shape in which the claws are absent. In particular, the subpixel averaging at points on the top of recognizer seems to depend on the presence or absence of claw feet at the same horizontal position on the opposite side of the structure, even though those are far away and shouldn't affect local considerations. I've been working on trying to pinpoint what is going on with this. @stevengj, as the designer of the sub-pixel averaging algorithm, do you have any guesses as to why subpixel averaging behaves differently for the recognizer and the simple inverted U? |
|
One of point-in-object, normal-to-object, or intersect-line-with-object must be giving a different result for the "U" and the "recognizer" shapes, up for a pixel at the top far from the claws where they really should be identical. So,
|
|
@HomerReid, you said something today about For example, the following test program, which computes the overlap of a box and a sphere, gives a continuously varying result as I vary the input argument #include <ctlgeom.h>
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char **argv)
{
geom_box b = {{0,-2,-2},{1,2,2}};
vector3 c = {0,0,0};
geometric_object o = make_sphere(NULL, c, 1);
double x = argc > 1 ? atof(argv[1]) : 0;
b.low.x += x;
b.high.x += x;
printf("overlap = %g\n", box_overlap_with_object(b, o, 1e-10, 1000000));
return 0;
} |
|
In libctl#20 I have overhauled the Meanwhile, images of the recognizer and pincer structures appear to be much improved. I'm continuing to investigate, but I think it's ready for more testing when you have a chance.
|
|
There is still a bug in the prism objects which produces different results when compared to an identical structure made out of blocks. This is demonstrated in the simulation script below which involves the "claw" structure from above created using either a prism ( Simulation Script import meep as mp
cell_size = mp.Vector3(5,5,5)
Si = mp.Medium(index=3.5)
h = 1.0
use_prism = True
if use_prism:
v = [mp.Vector3(-1.0,0,-1),
mp.Vector3(-1.0,0,-0.2),
mp.Vector3(-0.5,0,-0.2),
mp.Vector3(-0.5,0,-0.8),
mp.Vector3(1.0,0,-0.8),
mp.Vector3(1.0,0,0.8),
mp.Vector3(-0.5,0,0.8),
mp.Vector3(-0.5,0,0.2),
mp.Vector3(-1.0,0,0.2),
mp.Vector3(-1.0,0,1),
mp.Vector3(1.5,0,1),
mp.Vector3(1.5,0,-1)]
geometry = [mp.Prism(v, axis=mp.Vector3(0,1,0), height=h, material=Si)]
else:
geometry = [mp.Block(material=Si, size=mp.Vector3(2.5,h,2.0), center=mp.Vector3()),
mp.Block(material=mp.air, size=mp.Vector3(1.5,h+0.1,1.6), center=mp.Vector3()),
mp.Block(material=mp.air, size=mp.Vector3(0.5,h+0.1,0.4), center=mp.Vector3(-1.0,0,0))]
boundary_layers = [mp.PML(thickness=1.0)]
fcen = 1.0
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=0.2*fcen), component=mp.Ey, center=mp.Vector3())]
sim = mp.Simulation(resolution=20,
geometry=geometry,
cell_size=cell_size,
boundary_layers=boundary_layers,
sources=sources)
def print_stuff(sim):
p = sim.get_field_point(mp.Ex, mp.Vector3(-0.21,0.37,0.18))
print("ex:, {}, {}".format(sim.meep_time(), p.real))
sim.run(mp.at_beginning(mp.output_epsilon), mp.at_every(5,print_stuff), until_after_sources=0)Output (blocks) Geometry (blocks)
Output (prism) Geometry (prism)
|
|
Even with a convex structure and the subpixel averaging turned off, the results are significantly different when using either a prism ( Simulation Script import meep as mp
import math
cell_size = mp.Vector3(5,5,5)
Si = mp.Medium(index=3.5)
use_prism = True
if use_prism:
vertices = [mp.Vector3(-1,0),
mp.Vector3(-0.5,math.sqrt(3)/2),
mp.Vector3(+0.5,math.sqrt(3)/2),
mp.Vector3(+1,0),
mp.Vector3(+0.5,-math.sqrt(3)/2),
mp.Vector3(-0.5,-math.sqrt(3)/2)]
geometry = [mp.Prism(vertices, height=1.0, material=Si)]
else:
geometry = [mp.Block(material=Si, size=mp.Vector3(1,math.sqrt(3),1), center=mp.Vector3()),
mp.Block(material=Si, size=mp.Vector3(1,math.sqrt(3),1), center=mp.Vector3(),
e1=mp.Vector3(0.5,+math.sqrt(3)/2,0), e2=mp.Vector3(-math.sqrt(3)/2,0.5,0), e3=mp.Vector3(0,0,1)),
mp.Block(material=Si, size=mp.Vector3(1,math.sqrt(3),1), center=mp.Vector3(),
e1=mp.Vector3(0.5,-math.sqrt(3)/2,0), e2=mp.Vector3(+math.sqrt(3)/2,0.5,0), e3=mp.Vector3(0,0,1))]
boundary_layers = [mp.PML(thickness=1.0)]
fcen = 1.0
sources = [mp.Source(mp.GaussianSource(fcen, fwidth=0.2*fcen), component=mp.Ez, center=mp.Vector3())]
sim = mp.Simulation(resolution=20,
eps_averaging=False,
geometry=geometry,
cell_size=cell_size,
boundary_layers=boundary_layers,
sources=sources)
def print_stuff(sim):
p = sim.get_field_point(mp.Ex, mp.Vector3(-0.21,0.37,0.18))
print("ex:, {}, {}".format(sim.meep_time(), p.real))
sim.run(mp.at_beginning(mp.output_epsilon), mp.at_every(5,print_stuff), until_after_sources=0)Output (blocks) Geometry (blocks) Output (prism) Geometry (prism) The geometry looks similar for the two cases. The difference in the results may be due to the prism being shifted or stretched by a pixel or two relative to the identical structure created using blocks. This can be investigated by comparing the epsilon HDF5 file output for both runs (e.g., An isosurface plot of
|
|
After adjusting the position of the "claw"-shaped prism with the new Simulation Script import meep as mp
cell_size = mp.Vector3(5,5,5)
Si = mp.Medium(index=3.5)
h = 1.0
use_prism = True
if use_prism:
v = [mp.Vector3(-1.0,0,-1),
mp.Vector3(-1.0,0,-0.2),
mp.Vector3(-0.5,0,-0.2),
mp.Vector3(-0.5,0,-0.8),
mp.Vector3(1.0,0,-0.8),
mp.Vector3(1.0,0,0.8),
mp.Vector3(-0.5,0,0.8),
mp.Vector3(-0.5,0,0.2),
mp.Vector3(-1.0,0,0.2),
mp.Vector3(-1.0,0,1),
mp.Vector3(1.5,0,1),
mp.Vector3(1.5,0,-1)]
geometry = [mp.Prism(v, axis=mp.Vector3(0,1,0), center=mp.Vector3(), height=h, material=Si)]
else:
geometry = [mp.Block(material=Si, size=mp.Vector3(2.5,h,2.0), center=mp.Vector3()),
mp.Block(material=mp.air, size=mp.Vector3(1.5,h+0.1,1.6), center=mp.Vector3()),
mp.Block(material=mp.air, size=mp.Vector3(0.5,h+0.1,0.4), center=mp.Vector3(-1.0,0,0))]
sim = mp.Simulation(resolution=50,
eps_averaging=True,
geometry=geometry,
cell_size=cell_size)
sim.run(mp.at_beginning(mp.output_epsilon), until=0)We can generate images (shown below for two different viewing angles) of the difference in the permittivity values using Mayavi: Difference in Permittivity Values (with subpixel averaging)
The results look different when subpixel smoothing is turned off. Difference in Permittivity Values (without subpixel averaging)
|
|
Try it with |
|
When the prism has
|
|
With the blocks shifted as well, the difference in permittivity is:
|
|
Thanks for this useful test case. I believe the issues are resolved by libctl#23. I have created an updated version of the python code, which instantiates the TRON recognizer in two ways---as a single prism and as a union of Here is the revised python script, The script accepts three options:
Here's a shell script that runs the python code with both prism and block implementations and plots the two structures, plus their difference, in mayavi: As far as I can tell, the two structure plots are identical and the difference plot is all zeros. |
|
stevengj/libctl#23 does indeed produce a difference plot with all zeros when comparing the prism structure with a superposition of 5 import meep as mp
cell_size = mp.Vector3(5,5,5)
Si = mp.Medium(index=3.5)
h = 1.0
use_prism = False
if use_prism:
v = [mp.Vector3(-1.25,0,-1),
mp.Vector3(-1.25,0,-0.2),
mp.Vector3(-0.75,0,-0.2),
mp.Vector3(-0.75,0,-0.8),
mp.Vector3(+0.75,0,-0.8),
mp.Vector3(+0.75,0,+0.8),
mp.Vector3(-0.75,0,+0.8),
mp.Vector3(-0.75,0,+0.2),
mp.Vector3(-1.25,0,+0.2),
mp.Vector3(-1.25,0,+1),
mp.Vector3(+1.25,0,+1),
mp.Vector3(+1.25,0,-1)]
geometry = [mp.Prism(v, axis=mp.Vector3(0,1,0), center=mp.Vector3(), height=h, material=Si)]
else:
geometry = [mp.Block(material=Si, size=mp.Vector3(2.5,h,2.0), center=mp.Vector3()),
mp.Block(material=mp.air, size=mp.Vector3(1.5,h,1.6), center=mp.Vector3()),
mp.Block(material=mp.air, size=mp.Vector3(0.5,h,0.4), center=mp.Vector3(-1.0,0,0))]
sim = mp.Simulation(resolution=50,
geometry=geometry,
cell_size=cell_size,
eps_averaging=False)
sim.run(mp.at_beginning(mp.output_epsilon), until=0)
Difference Plot of Prism vs. Blocks
|
|
Just to clarify, the issue being reported here concerns a difference between two different ways of representing the same structure with the legacy I noted in the comments on libctl#23 that the legacy implementation does not appear to be meticulous in adjudicating questions regarding points lying on object boundaries, and perhaps what you are observing is a manifestation of that. In either case, the question is only relevant for cases in which grid points happen to lie on object boundaries---a matter of coincidence that is purely an artifact of specific discretization and should vanish (even at the same resolution!) upon slight (subpixel) displacements of the objects to shift boundaries off of grid points. My interpretation of the bedrock philosophy of meep is that the code guarantees convergence to correct results with increasing resolution, but doesn't guarantee equality of results obtained on two shifted grids of the same resolution. Do we want to revise that policy to offer new, stronger, promises regarding effects like the ones you've observed? Note that the inconsistency you've observed is in the legacy Either way, since the original issue seems to have been resolved and the discussion seems to have morphed into larger questions, I propose closing the issue, though I'd like to give the original poster a chance to determine whether the issue has been resolved to his satisfaction. |
|
Thanks for the libctl update Homer! I've rerun my test suite using libctl#23, and can confirm that the volume errors on the "recognizer" shapes are no longer present. I also ran the other test case above (the strip to slot waveguide converter) and it looks good and without the previous errors, so I agree that closing it would be appropriate. I do notice some strange discretization effects (in terms of surface pixels) in other examples with multiple prisms coincident at the top/bottom surfaces, but I think perhaps this is out of the scope of the current bug report. Thanks again! |
|
Let's close this since the original problem is gone. If additional problems appear, let's open new issues for those. |
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.
The text was updated successfully, but these errors were encountered: