VertexOnlyMesh and Extruded Meshes

[ReubenHill]

So I’m trying to get VertexOnlyMesh to play nicely with extruded meshes. One issue is getting the runtime quadrature to play nicely. I can get VertexOnlyMeshTopology.cell_parent_cell_list to produce the correct output for UnitSquareMesh(3, 3, quadrilateral=True) and ExtrudedMesh(UnitIntervalMesh(3), 3) - in each case the numbering varies but the below works:

def test_extruded_cell_parent_cell_list():
    # If we make a function space that has 1 dof per cell, then we can use the
    # cell_parent_list directly to see if we get expected values. This is a
    # carbon copy of tests/regression/test_locate_cell.py

    ms = UnitSquareMesh(3, 3, quadrilateral=True)
    mx = ExtrudedMesh(UnitIntervalMesh(3), 3)

    # coords at locations from tests/regression/test_locate_cell.py - note that
    # we are not at the cell midpoints
    coords = np.array([[0.2, 0.1], [0.5, 0.2], [0.7, 0.1], [0.2, 0.4], [0.4, 0.4], [0.8, 0.5], [0.1, 0.7], [0.5, 0.9], [0.9, 0.8]])

    vms = VertexOnlyMesh(ms, coords)
    vmx = VertexOnlyMesh(mx, coords)
    assert vms.num_cells() == len(coords)
    assert vmx.num_cells() == len(coords)
    assert np.equal(vms.coordinates.dat.data_ro, coords).all()
    assert np.equal(vmx.coordinates.dat.data_ro, coords).all()

    # set up test as in tests/regression/test_locate_cell.py - DG0 has 1 dof
    # per cell which is the expression evaluated at the cell midpoint.
    Vs = FunctionSpace(ms, 'DG', 0)
    Vx = FunctionSpace(mx, 'DG', 0)
    fs = Function(Vs)
    fx = Function(Vx)
    xs = SpatialCoordinate(ms)
    xx = SpatialCoordinate(mx)
    fs.interpolate(3*xs[0] + 9*xs[1] - 1)
    fx.interpolate(3*xx[0] + 9*xx[1] - 1)

    # expected values at coordinates from tests/regression/test_locate_cell.py
    expected = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])
    assert np.allclose(fs.at(coords), expected)
    assert np.allclose(fx.at(coords), expected)
    assert np.allclose(fs.dat.data[vms.cell_parent_cell_list], expected)
    assert np.allclose(fx.dat.data[vmx.cell_parent_cell_list], expected)

The above numbering is the same as you get if you ask the extruded mesh to do locate_cell on a given point.

I use VertexOnlyMesh.cell_parent_cell_list to make a pyop2 Map to the parent mesh cell_set - of course the cell_set for an extruded mesh is a pyop2 ExtrudedSet: this uses the base mesh cell_set and the given extruded layers so I don't know if I should expect the locate_cell list of numbers to play nicely with this.

When I interpolate, I make a pyop2 ComposedMap from the VertexOnlyMesh.cell_parent_cell_list to FunctionSpace.cell_node_map which I then use in my par loop. @ksagiyam or @connorjward , can you advise at all?

[connorjward]

So my understanding of what you are doing is you want to have a parloop that does something like:

for pt in vertex_only_mesh:
    // pack the data found in the closure (e.g.) of the cell which contains pt
    for a in arity:
        t0[a] = extr_dat[pt_to_dofs[pt*arity + a]]

    // compute
    clever_kernel(&t0)

    // unpack

I think what you should look to do is to compute and use a new Map where you:

1. Find the right cell number with locate_cell (should work fine I think)
2. Identify the right base cell and column height that correspond to this cell
3. Store in the new map something like:

for pt in vertex_only_mesh:
    base_cell, height = unflatten_locate_cell(pt)

    for a in arity:
        new_map[pt*arity + a] = extr_map[base_cell*arity+a] + extr_map.offsets[a]*height

Basically I think things get easier if you stop thinking about the Dat being extruded and just address it using standard maps. Otherwise I think you'd need a strange type of map that mapped from a vertex-only mesh cell to a 2-tuple of (base_cell, height).

I could be totally misunderstanding what you want, and I could be wrong, but I hope this helps.

[ReubenHill]

It works! Thanks for the help. Draft PR here, hopefully it passes all tests #2750 🤞