single-cell periodic quad work

[ThighamOxford]

Support single-cell doubly periodic quadrilateral closure

Handle the four-entity DMPlex closure that occurs for a single quadrilateral cell periodic in both directions.

DMPlex returns the unique cell, two edges, and one vertex. Expand these to the nine local entities expected by Firedrake by repeating each edge twice and the vertex four times, analogous to the existing nclosure == 6 handling.

Description

This is a draft implementation adding closure ordering support for a single quadrilateral cell that is periodic in both directions.

DMPlexGetTransitiveClosure returns only the four unique entities in this case: one cell, two edges and one vertex. Firedrake's quadrilateral closure ordering expects the nine local entities of the reference quadrilateral.

This change handles the nclosure == 4 case analogously to the existing nclosure == 6 handling, repeating the two edges and single vertex to construct the expected nine-entry local closure.

With this change, PeriodicRectangleMesh(1, 1, ..., direction="both", quadrilateral=True) successfully constructs and produces the expected cell closure.

Still to do / discussion

The continuous coordinate field is degenerate for a one-cell periodic direction. I have verified locally that constructing DQ1 coordinates and passing them through _postprocess_periodic_mesh gives the correct localized geometry. My test problem is a 2π×2π one-cell torus on which I integrate 1 to get an area of 4π.
Without the DQ coordinate postprocessing, I think the 1 × 1 doubly-periodic PeriodicRectangleMesh has only a single coordinate value (0, 0), since all four corners are identified with the same topological vertex. This results in a geometrically collapsed cell (zero area). I think this also holds for the 1 x ny, and nx x 1 doubly periodic cases where ny, nx >1.

I'm opening this as a draft while working out the preferred way to integrate the localized-coordinate handling into PeriodicRectangleMesh, and would welcome guidance on whether reusing _postprocess_periodic_mesh is the intended approach.

Equally I might be incorrect!

[connorjward]

I don't think you need to do anything around DG coordinates. It should just work.

Could you add a test for this? Probably in https://github.com/firedrakeproject/firedrake/blob/main/tests/firedrake/regression/test_mesh_generation.py using existing tests as inspiration.

[ThighamOxford]

I added a test for this problem - when I run it myself it fails (returns an integral of zero). I still think there is an issue with the coordinates of the nodes.

[connorjward]

I added a test for this problem - when I run it myself it fails (returns an integral of zero). I still think there is an issue with the coordinates of the nodes.

Thanks. I can investigate this myself now.

Before that, could you push your branch to the main firedrakeproject repository instead of your fork? That way I can edit things.

[ThighamOxford]

Thanks Connor - I've done that now

[connorjward]

I have pushed a bunch of changes (you need to change the branch used in this PR to the one that lives in this repository). I got it to work for the case you originally had but am now seeing failures for some more general cases (quadrilateral=False and direction!="both"). I've also had to put in some parallel fixes (which may solve #5371?).

The crux of the issue is that PETSc does not handle periodic coordinates well. The main algorithm for determining DG coordinates fails when there is only a single cell. I therefore have to do things by hand - not hard, just annoying.

I've kind of run out of steam but I can revisit later this week. I think the remaining things to do are:

• Open a PETSc issue explaining the problem at their end (I should do this).
• Either fix the issues with the other refinement types, or: raise NotImplementedError, xfail the tests, and create an issue.
• Make sure the code is all cleaned up and documented