What is the point of RiemannianMetric._set_plex_coordinates()

[stephankramer]

I'm trying to understand the logic of what we do in _set_plex_coordinates(), but if my understanding is correct (and that's what I'd like to check here really) the reordering we do there is superfluous. The reason for my question, is that it needs fixing for higher order (curved) meshes - but afaics I could also "fix it" by just deleting it.

So I think there are 3 numberings at play here:

1. The ordering of the vertices as "topological points" which are contiguously numbered (as a "stratum") from vStart to vEnd
2. The ordering in which the coordinates appear in the vector that you get from getCoordinatesLocal, and there's mapping from 1. to 2. provided by the section of the coordinate dm. So this is an ordering used by petsc internally, which in principal might be different from 1
3. The ordering used by Firedrake when storing the DOF of P1 discretisations in a vector. Again the mapping from 1. to 3. is provided by a section, which in _set_plex_coordinates() we (re)create by calling mesh.create_section() asking for the appropriate ndofs at each vertex. This ordering determines the order of e.g. mesh.coordinates.dat.data[:] (assuming we have a linear, not a curved mesh) and the ordering of the .dat.data[:] of the metric itself (just with a different ndofs per vertex).

So what it appears is that _set_plex_coordinates() changes the ordering of 2. to be the same as 3.

Now if you look in the petsc code of src/dm/impls/plex/adaptors and src/dm/impls/plex/plexmetric.c the only place where the local coordinate vector is actually used, is in the interface to (par)mmg where it's also using the coordinate dm section, and then if you look at src/dm/impls/plex/adaptors/mmg/mmgadapt.c:

  for (v = 0; v < vEnd - vStart; ++v) {
    PetscCall(PetscSectionGetOffset(coordSection, v + vStart, &off));
    for (i = 0; i < dim; ++i) vertices[dim * v + i] = PetscRealPart(coords[off + i]);
  }

so the resulting vertices array that's passed to mmg is always in the same ordering as 1. (i.e. the naive contiguous vertex ordering) - and that has to be the case because when we pass a cell with vertices [0, 1, 2] to mmg (and these numbers are straight topological point numbers) then mmg will expect those vertices to correspond the 0th, 1st and 2nd entries in the vertices array. So that means the interfaces doesn't care what numbering we use for 2. The code in the parmmg interface is a bit more complicated, as it makes a selection of vertices to pass on to parmmg, but still the same things applies: we can use any ordering of the coordinates vector (ordering 2) as long as the coordinate dm section is consistent with it.

If you look at how the metric is passed on however, it doesn't appear to use any section at all! It seems to completely assume that the ordering in the metric vector is the same as the topological point numbering of the vertices (i.e. ordering 1) - and so that's why I think we need to use this to_petsc_local_numbering(): it reorders the metric vector (it's .dat.data[:]) which will be in ordering 3. into ordering 1.

But all of that seems to imply that there really is no point of trying to make ordering 2. (the order in the local coordinate vec) to be the same as that in 3. (Firedrake's P1 ordering) - because petsc never sees that ordering 3. when we go into (par)mmg, and the interface should work regardless of what ordering 2. we use.

So the only reason/logic for _set_plex_coordinates() I can see is that maybe it is necessary if you want to use petsc's normalisation routine DMPlexMetricNormalize() ? - as that involves a FEM integral which does also use the coordinate vector - but of course we have our own implementation of that normalisation, where all concerns about ordering are Firedrake's. Any of the DMPlexMetric functions that we do actually use (determinants,intersection, etc.) does not actually use the petsc coordinate vec AFAICS

Sorry, this is a bit of blurb, just organizing my thoughts as well by writing it out...