Bug in `cellmask!` after `uniform_refine` a 3D grid

[pjaap]

MWE: create a cube and edit the region of the upper half

using ExtendableGrids, GridVisualize, GLMakie

The following code

grid = uniform_refine(grid_unitcube(Tetrahedron3D))
cellmask!(grid, [0.0, 0.0, 0.5], [1.0, 1.0, 1.0], 2) # upper half should be = 2
gridplot(grid, Plotter = GLMakie, xplanes = [0.0001])

results in

...which is wrong (see the missing green tetrahedron)
While creating the grid "directly"

grid = simplexgrid(0:0.5:1, 0:0.5:1, 0:0.5:1)
cellmask!(grid, [0.0, 0.0, 0.5], [1.0, 1.0, 1.0], 2) # upper half should be = 2
gridplot(grid, Plotter = GLMakie, xplanes = [0.0001])

results in

...which looks correct

@hadjimy fyi (this causes the weird demo results in our session)

[pjaap]

to add: this is not a round-off problem. Using cellmask!(grid, [0.0, 0.0, 0.49], [1.0, 1.0, 1.0], 2) yields the same wrong result.

something weird is going on in the CellNodes of the refined grid:
grid[Coordinates][:, grid[CellNodes][:, 22]]:

3×4 Matrix{Float64}:
 0.0  1.0  0.5  0.5
 0.0  0.5  0.5  0.5
 0.5  1.0  0.5  1.0

hence, cell 22 has x-coordinates 0 and 1. this cannot happen in a refined grid, where each tetrahedron has size 0.5.

[pjaap]

Oh dear.

@Liameloy and I leaned a lot today. This is not a bug. I simply had a misconception about refinements of a tetrahedron.

For future readers, and future me: A tetrahedron cannot be split into 8 congruent tetrahedra. At the corners, we have 4 congruent ones with half length, but the center remainder in formed by an octahedron. In ExtendableGrids, a refinement into 4 tretrahedra is chosen, such that one common edge is used. And this edge can have length 1 in the picture above.

By using yplanes and zplanes above, the other larger tetrahedra are being revealed.

@j-fu Potentially, uniform refinements in this fashion may lead to poorly shaped cells. Was there a discussion to split the central octahedron into 8 congruent tetrahedra by adding a new center node?

[j-fu]

No, the refinement code has been written by @chmerdon .

[chmerdon]

Yes, I did this without much discussion. Just wanted to have something that was working fast at that time ^^;;

I once started a pull request to unify the split_into and uniform refinements into one function and introduce refinement types to steer the refinement. But it somehow never got finished. Maybe one could start a new attempt and also integrate alternative refinement that preserve the mesh quality.

[chmerdon]

Does it really decrease the mesh quality, though? It maybe depends on which of the two diagonals in the central octagon is chosen for the refinement. Currently no smart choice is implemented.