How to ensure a line is included in the mesh? Possible source of RuntimeError: Factor is exactly singular

[duarte-jfs]

I am trying to include a LineString with the purpose of using it for a line integral. However, that line must go through other polygons, and I am not understanding when it is included and when it is not. When it appears to be included, the modes can be computed. When it is not, then it raises an error of RuntimeError: Factor is exactly singular. Any clues how to solve this?

Code

h = 0.2 #CHANGE HERE
W = 5

box1 = box(0,1,1,1+h)
box2 = box(2,0,3,1)
box3 = box(0,0,2,1.)
background = box(-W+1,-W+0.5,W+1,W+0.5)

line = LineString([(0,(box1.bounds[1]+box1.bounds[3])/2),
                  (1,(box1.bounds[1]+box1.bounds[3])/2),
                  (2,0.5),
                  (3, 0.5)])

polygons = OrderedDict(box1 = box1, 
                       box2 = box2,
                       line = line,
                       box3 = box3,
                       background = background
                       )

resolutions = dict(
    box1 = {'resolution' : 1, 'distance': 1},
    box2 = {'resolution' : 1, 'distance': 1},
    box3 = {'resolution' : 0.2, 'distance': 1},
    line = {'resolution' : 0.1, 'distance': 0.1},
    background = {'resolution': 1, 'distance':1}
)



mesh = from_meshio(mesh_from_OrderedDict(polygons, 
                                         resolutions))

basis = Basis(mesh, ElementTriP0(), intorder=4)

epsilon = basis.ones(dtype = complex)

        
elements = basis.get_dofs(elements = 'box3') 
epsilon[elements] = 3**2
        
fig = mesh.draw()
fig.set_axis_on()


ax = plt.gca()

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), alpha = 0.5)
    else:
        ax.plot(*polygon.xy, alpha = 0.5)
        
modes = compute_modes(
    basis,
    epsilon,
    1550,
    num_modes = 1
)

modes[0].show('E', colorbar = True)
# fig = plt.figure()
# ax = fig.add_subplot(111)

Running the code above with ``h=0.2` gives the error mentioned above, and when looking at the mesh superimposed with the polygons boundaries we see that the line is not included (green line):

If we now change h=1 then the mode is computed just fine and we now see that the line is clearly included in the mesh:

I'm unsure how to interpret this

EDIT

I have narrowed down the error to a precision one. The problem stems from the function break_line_ in mesh.geometry as it tries to split a line my intersection when intersection.geo_type = 'Point' via:

line = linemerge(split(line, intersection))

This works in some cases and not in others. Yet the intersection is always found. For now I have fixed it through:

line = snap(line, intersection))

I'm unsure if this is the best way to do it, but it caused a significant problem in my case as I was interested in calculating a line integral but since half the line was not being meshed, the calculated integral was far lower than expected.

[HelgeGehring]

Did you try to put the line at highest mesh order?

polygons = OrderedDict(
                       line = line,
                       box1 = box1, 
                       box2 = box2,
                       box3 = box3,
                       background = background
                       )

Lines should have the highest priority to not be overwritten (the boxes will still be there anyway)

But even without, the line should be there outside of the two boxes with higher priority 🤔

@simbilod what do you think? Is this the fix or is there a problem with the linemerge?

[duarte-jfs]

Yes, I tried that order as well, but the error just kept on coming back :/

[HelgeGehring]

Hmm, we should fix that :/

If you use snap, does it always work for you? I'd expect that it tries to snap the points of the line to other points on the grid?

@simbilod will have more insight in the meshing, but he probably won't have too much time the next days

[duarte-jfs]

So far, I believe it has, but I haven't checked thoroughly. I'm using it in a very big loop (1000 iterations) and a few (probably less than 20) have given an error, which stems from the presence of the the line that I'm trying to mesh. But I'm not sure if it is related to the the exact same problem as above. If I find out I will update here.

[HelgeGehring]

Sounds good!