Bug: Periodic cartesian meshes

[JordiManyer]

So there is a bug with periodic cartesian models in GridapP4est, where the periodicity is not correctly applied to the p4est mesh. See the following MWE:

using Gridap, Gridap.Geometry
using PartitionedArrays, GridapDistributed
using GridapP4est

ranks = with_mpi() do distribute
  distribute(LinearIndices((1,)))
end

cmodel = CartesianDiscreteModel((0,1,0,1,0,3),(1,1,3);isperiodic=(false,false,true))

model = OctreeDistributedDiscreteModel(ranks,cmodel,0)
m = PartitionedArrays.getany(local_views(model))

# These two are different
ctopo = get_grid_topology(cmodel)
Geometry.get_faces(ctopo,3,0) # Periodic, with 12 nodes
Geometry.get_vertex_coordinates(ctopo)

topo = get_grid_topology(m) 
Geometry.get_faces(topo,3,0) # Not periodic, with 16 nodes

The problem comes from the fact that the p4est connectivity is computed from the model cell nodes ids and node coordinates (see here). However, in a CartesianDiscreteModel those seem to NOT have periodicity. I.e reshape(collect(get_cell_node_ids(cmodel)),(3,)) !== Geometry.get_faces(ctopo,3,0) for some reason.

The periodic cell-to-node map can only be obtained from the topology (see fix). I don't know why this happens, and it might also be a bug (?).

---

EDIT: So the issue is that the model cannot be periodic, only the topology can. Otherwise it breaks the cell-maps, which are FE functions based on the node coordinates. For reference, with the current changes we have that this code gives the wrong answer:

φ3 = get_cell_map(m)[3]
pt = VectorValue(0.5,0.5,0.5)
φ3(pt) # VectorValue{3, Float64}(0.5, 0.5, 1.0) WRONG!

This means the fix is a lot more involved, since the connectivity for the grid and the topology are different (I believe). So my fix does not work.

[JordiManyer]

Another MWE, now in 2D, illustrating my latter point:

m1 = CartesianDiscreteModel((0,3,0,3),(3,3);isperiodic=(true,true))
m2 = CartesianDiscreteModel((0,3,0,3),(3,3);isperiodic=(false,false))

t1 = get_grid_topology(m1)
t2 = get_grid_topology(m2)

# These two are very different
Geometry.get_faces(t1,2,0)
Geometry.get_faces(t2,2,0)

# These two are the same
get_cell_node_ids(m1)
get_cell_node_ids(m2)

The cell-to-node map are completely different between the grid and the topology...

[JordiManyer]

Update: I now believe that P4est cannot correctly handle periodic meshes on it's own...

The current fix sets the p4est connectivity to correctly reflect the periodicity. This seems to cause p4est to incorrectly compute the new vertex coordinates in the cells that contain periodic nodes. It seems to simply use a linear combination of the coarse cell vertex coordinates to determine the new vertex coordinates (which is wrong in cells containing periodic nodes).

[amartinhuertas]

Update: after talking in Slack, it seems to be possible. in p4est_connectivity we have vertices (geometry) and corners (topology). At present, from GridapP4est, we are treating these two indistinguishably.

[JordiManyer]

Solved in #61 , closing issue.