StructuredTriangleGrid strange behaviour

[alessiofumagalli]

Dear all, I found a strange behaviour related to the function StructuredTriangleGrid. When N <= 15 is set the grid looks fine, but from 16 the ordering of the elements are not coherent with what we are expecting. Also, this has some implications with the RT0 discretization. I leave here a minimal example:

N = 15 # 16
g = pp.StructuredTriangleGrid([16]*2)
g.compute_geometry()

pp.plot_grid(
    g,
    plot_2d=True,
    info="fn",
    alpha = 0,
    fig_size = (30, 24)
)

Thanks a lot for your kind help :)

[keileg]

Hi Alessio. Can you be more specific on what is unexpected to you? Also, the plot you generate show only face and node indices, from what I can tell.

[alessiofumagalli]

For N <= 15 the cells, faces and nodes are numbered by row from the bottom of the domain. However, starting from 16 this ordering is not kept anymore. I feel something went wrong. See for example near the top left corner of the following grid (to be zoomed)

[keileg]

You are right, something has changed here. I traced the issue down to pp.utils.setmembership.unique_rows. That function uses a strange trick that I found on stackexchange back in the days. Presumably, something in numpy has changed which distorts the data arrays when uniquifying. To be clear, the function still functions as promised, it just behaves in a different way than one could expect. Thus code that depends on the ordering is applying unwarranted tacit assumptions.

Since this function was written, numpy has added an axis option to unique, and thus I believe the opaque code can be changed to (roughly)

def unique_rows(data: np.ndarray):
    return np.unique(data, axis=0, return_index=True, return_inverse=True)

At the first glance, this gives a more natural ordering of the returned values, thus of the face-node relation. I tried to introduce this, but got a lot of error messages. Incidentally, the majority of the errors were from tests that you have written...

I think that the long term solution is to change the implementation of unique_rows and make the needed changes to the rest of the code. In the shorter term (at least two weeks), I do not have the capacity to look further into this (side note: we are doing a major refactoring of the test suite next week, as detailed in newly opened issues; you may in particular want to have a look at #967, both list of tests and comment below). So, I will leave this issue open for now, and return to it in due course.