Skip to content

[BUG] utils/geom.c: unecessary error raise in init_prism when two consecutive vertices are alligned with the centroid. #61

@Niccolo-Marcucci

Description

@Niccolo-Marcucci

When initiating a prism in meep i was getting the following error:

-----------
Initializing structure...
non-coplanar vertices in init_prism

even though all the prism vertices were such that vertices[i].z = 0.0. I guess this is the error raised at line 2598 in the file utils/geom.c.
I then noticed that, for checking the prism plane, init_prism() checks the normals of all the triangles composed of two consecutive vertices and the centroid of the prism face.

Could it be that having consecutive vertices aligned with the centroid would give rise to an error exactly because the resulting triangle would have null area? (thus the normal might be affected by numerical error).

Steps to reproduce the bug:

import meep as mp
import numpy as np
import matplotlib.pyplot as plt
sim = mp.Simulation(
            cell_size = mp.Vector3(6,6,6),
            geometry = [],
            sources = [],
            resolution = 20,
            boundary_layers = [],
            dimensions = 3,
            symmetries = [],
            force_complex_fields = False,
            eps_averaging = False)

# frame-like prism
vx = [1,1,-1,-1,1, 1, 2,2,-2,-2,2,2]
vy = [0,1,1,-1,-1, 0, 0,-2,-2,2,2,0]
v = [ mp.Vector3(vx[i],vy[i],0) for i  in range(len(vx))]
centroid =  sum(v, mp.Vector3(0)) * (1.0 / len(v))

c1 = mp.Prism(vertices = v,
              height = .3,
              axis = mp.Vector3(0,0,1),
              center = mp.Vector3() + centroid,
              material=mp.Medium(epsilon=2))
sim.geometry.append(c1)

sim.init_sim()

simsize = sim.cell_size

fig = plt.figure(dpi=200)
plot = sim.plot2D( output_plane=mp.Volume(size=mp.Vector3(simsize.x,simsize.y)),
                labels=True,
                eps_parameters={"interpolation":'none',"cmap":'gnuplot', "vmin":'0'})

Things that make the error message disappear

  • change the y-array to vy = [0,1,1,-1,-1, 0, 0.001,-2,-2,2,2,0]; in this way there are no consecutive points aligned with the centroid
  • comment the line relative to the center assignment in the Prism() block.
  • weirdly so, some velues of thickness make the error disappear. For instance:
    -- 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 0.9, 0.03, 0.02, 0.01, ... in general, decimals do not work
    -- 0.7 works
    -- any number in the form 2**N or 1+2**N, with N negative integer number, works.
    -- integer number work

Although my explanation seem reasonable, i really do not understand why changing the thickness should affect the raising of the error.
I'm unable to explore further.
Would it be possible to solve the problem by checking the area of the triangle, or checking for aligned vertices before computing the normal?

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions