Meshing error: "Exception: Start and end points of line should be different"

[duarte-jfs]

I am unable to understand the source of the error below:

# box1 = box(0, 0, 1, 1)
# box2 = box(1,-2,5,5)
# box3 = box(0,1,1,2)

box1 = box(-1.0, 1.0, 1.0, 1.454) 
box2 = box(1.0, 0.2, 17.0, 2.6540000000000012) 
box3 = box(-1.0, 1.454000000000001, 1.0, 1.5080000000000011) 

surface = box1.buffer(0.5).exterior
polygons = OrderedDict(box1 = box1,
               box2 = box2,
                box3 = box3)

resolutions = OrderedDict(box1 = {'resolution': 0.2, 'distance':1},
                   box2 =  {'resolution': 0.5, 'distance':1},
                   box3 =  {'resolution': 0.8, 'distance':1})

fig = plt.figure()
ax = fig.add_subplot(111)

for key, polygon in polygons.items():
    if type(polygon) not in [LineString, MultiLineString]:
        x,y = polygon.exterior.xy
        ax.plot(np.asarray(x),np.asarray(y))

    else:
        if type(polygon)==MultiLineString:
            for linestr in polygon.geoms:
                x,y = linestr.xy
                ax.plot(np.asarray(x),np.asarray(y), color = 'red')
        else:

            x,y = polygon.xy
            ax.plot(np.asarray(x),np.asarray(y), color = 'red')
            
mesh = from_meshio(mesh_from_Dict(polygons, resolutions))


mesh.draw()

Which gives:

Exception                                 Traceback (most recent call last)
Cell In[513], line 36
     33             x,y = polygon.xy
     34             ax.plot(np.asarray(x),np.asarray(y), color = 'red')
---> 36 mesh = from_meshio(mesh_from_Dict(polygons, resolutions))
     39 mesh.draw()

File ~\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\femwell\mesh\mesh.py:85, in mesh_from_Dict(shapes_dict, resolutions, default_resolution_min, default_resolution_max, filename, gmsh_algorithm, global_quad, verbose)
     83 meshtracker = MeshTracker(model=model)
     84 for i, polygon in enumerate(shapes_tiled_list):
---> 85     meshtracker.add_xy_surface(polygon, i, physical=False)
     86 # Tag physicals
     87 # TODO integrate in meshtracker
     88 surface_name_mapping = {}

File ~\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\femwell\mesh\meshtracker.py:188, in MeshTracker.add_xy_surface(self, shapely_xy_polygons, label, physical)
    185 for vertex in shapely.geometry.MultiPoint(shapely_xy_polygon.exterior.coords).geoms:
    186     # gmsh_point = self.add_get_point(vertex, label)
    187     exterior_vertices.append(vertex)
--> 188 channel_loop = self.xy_channel_loop_from_vertices(exterior_vertices, label)
    190 # Create and log surface
    191 gmsh_surface = self.model.add_plane_surface(channel_loop, holes=hole_loops)

File ~\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\femwell\mesh\meshtracker.py:77, in MeshTracker.xy_channel_loop_from_vertices(self, vertices, label)
     75 edges = []
     76 for vertex1, vertex2 in [(vertices[i], vertices[i + 1]) for i in range(len(vertices) - 1)]:
---> 77     gmsh_line, orientation = self.add_get_xy_segment(vertex1, vertex2, label)
     78     if orientation:
     79         edges.append(gmsh_line)

File ~\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\femwell\mesh\meshtracker.py:121, in MeshTracker.add_get_xy_segment(self, shapely_xy_point1, shapely_xy_point2, label)
    119     self.xy_segments_secondary_labels[index] = label
    120 else:
--> 121     gmsh_segment = self.model.add_line(
    122         self.add_get_point(shapely_xy_point1),
    123         self.add_get_point(shapely_xy_point2),
    124     )
    125     self.shapely_xy_segments.append(
    126         shapely.geometry.LineString([shapely_xy_point1, shapely_xy_point2])
    127     )
    128     self.gmsh_xy_segments.append(gmsh_segment)

File ~\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\pygmsh\common\geometry.py:77, in CommonGeometry.add_line(self, *args, **kwargs)
     76 def add_line(self, *args, **kwargs):
---> 77     return Line(self.env, *args, **kwargs)

File ~\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\pygmsh\common\line.py:25, in Line.__init__(self, env, p0, p1)
     23 assert isinstance(p0, Point)
     24 assert isinstance(p1, Point)
---> 25 id0 = env.addLine(p0._id, p1._id)
     26 self.dim_tag = (1, id0)
     27 self.dim_tags = [self.dim_tag]

File ~\AppData\Local\anaconda3\envs\phd_1\lib\site-packages\gmsh.py:6792, in model.occ.addLine(startTag, endTag, tag)
   6786 api_result_ = lib.gmshModelOccAddLine(
   6787     c_int(startTag),
   6788     c_int(endTag),
   6789     c_int(tag),
   6790     byref(ierr))
   6791 if ierr.value != 0:
-> 6792     raise Exception(logger.getLastError())
   6793 return api_result_

Exception: Start and end points of line should be different

Yet, when uncommenting the upper 3 lines, which is an identical geometry as before the code runs without any issue. What am I missing?

EDIT

Right after posting the question I narrowed down the problem, but I still don't know how to solve it nicely. It appears to be a numerical error. When switching between the two sets of lines of code, the top one runs just fine, while the bottom one doesn't. The only ifference is the numerical precision of one value. How can we avoid this error?

box1 = box(-1, 1, 1, 1.454)
box2 = box(1,0.2,17,2.654)
box3 = box(-1,1.454,1,1.508)

# box1 = box(-1, 1, 1, 1.454)
# box2 = box(1,0.2,17,2.654)
# box3 = box(-1,1.4540000000001,1,1.508)

[HelgeGehring]

@simbilod do you have an idea?

[duarte-jfs]

I found this solution: https://stackoverflow.com/questions/28028910/how-to-deal-with-rounding-errors-in-shapely

And have used the lines:

for key in polygons.keys():
    polygons[key] = shapely.set_precision(polygons[key], 1e-4)

so far so good, but not sure if it'll give more problems down the line

[simbilod]

Right, gmsh has its own point precision

In the meshwell / gdsfactory updated versions of these, I do internally round the coordinates of the input polygons (usually to dbu precision) as you show to avoid this