NanoComp / libctl

Guile-based library implementing flexible control files for scientific simulations
GNU General Public License v2.0
18 stars 23 forks source link

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

Open Niccolo-Marcucci opened 2 years ago

Niccolo-Marcucci commented 2 years ago

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

Yaraslaut commented 1 year ago

In geometry you have following setup: image Where point 1 is centroid (automaticaly calculated) and then points 2 and 3, during prism initialization normal vector of the triangle that a created from two vectors (1-2) and (1-3) highly sensible to numerical innacurecy of arithmetic operations and strictly speeking this two vector are parallel, this is the reason why such prism is not initilized properly and has weird beheviour while changing some numbers.