remove allocations from cell finder

[chmerdon]

The cell finder causes allocations, which affects performance of PointEvaluator or lazy_interpolate in ExtendableFEM.
So far I fixed one line that was causing 2 allocations in each call, another allocation seems to be hidden in (the call of) update_trafo! that I could not fix.

[j-fu]

If one writes

    L2G::L2GTransformer{Float64, Int32, Tetrahedron3D, Cartesian3D} = L2GTransformer(xgrid[CellGeometries][cell], xgrid, ON_CELLS)

the allocation goes away.

"Informally, a function is type stable if the type of the output depends only on
the types of the inputs, not their values"
There is now a whole PhD thesis on type stability in Julia: https://janvitek.org/pubs/artem-phd23.pdf .

I got this from Gemini after an hour staring at the code... The constructor of L2GTransformer is not type stable. I have to admit that this was not fully clear to me, though it is logical. Essentially, we might try to pass the CoordinateSystem along with the ElementGraph and see if this fixes it.

IMHO getting type stability solved would be worth a breaking change. But my gut feeling is that there are more of these cases.

[Da-Be-Ru]

If one writes

    L2G::L2GTransformer{Float64, Int32, Tetrahedron3D, Cartesian3D} = L2GTransformer(xgrid[CellGeometries][cell], xgrid, ON_CELLS)

the allocation goes away.

I cannot confirm this. We could also not confirm @chmerdon's allocation in update_trafo! on my machine.
But the changes so far do eliminate at least 2 allocations on my end.

Since a downstream PR is waiting on this to be merged, I'd suggest doing that and opening two separate issues on the L2GTransformer constructor and the remaining allocation hunt.

So can we merge this, @j-fu?

[j-fu]

I would agree. But if there is time: I thought a bit more on this. Adding a type stable method

function L2GTransformer(
        EG::Union{Type{<:Tetrahedron3D}, Type{<:Parallelepiped3D}},
        grid::ExtendableGrid{Tv, Ti},
        AT::Type{<:AssemblyType},
        CT::Type{<:AbstractCoordinateSystem}
    ) where {Tv, Ti}
    A = zeros(Tv, size(grid[Coordinates], 1), dim_element(EG))
    b = zeros(Tv, size(grid[Coordinates], 1))
    return L2GTransformer{Tv, Ti, EG, CT}(
        0,
        false,
        grid[Coordinates],
        grid[GridComponentNodes4AssemblyType(AT)],
        grid[GridComponentVolumes4AssemblyType(AT)],
        A,
        b,
        zeros(Tv, 3, 3),
        0
    )
end

could help to solve this. However one would have to do more at the call site as well.

[chmerdon]

I don't understand this yet. You suggested a new constructor that, as far as I can see, produces the same object as before, so where is the difference? I think, the allocations happen with update_trafo! (or not at all, since @Da-Be-Ru does not see them). The L2GTransformers used here for each geometry were already constructed only once and stored in the CellFinder before gFindLocal is called. (The problem in gFindLocal is that the geometry might change while traversing the grid, so one cannot use a single geometry and might have to switch to another L2GTransformer in the while loop. But maybe one can write a version for grids that only have a single unique geometry, so that one can use a barrier function approach...)

Anyhow, I am in favor of merging this PR now and open a new one for the other ideas.