.at() for periodic extruded meshes

[maxeruth]

Hi all, I'm just starting to learn Firedrake. I am interested in solving a PDE on the solid torus, i.e. $\Omega = D \times S^1$ where $D$ is the unit disk and $S^1$ is the circle. From my understanding, the natural method for doing this in Firedrake appears to be a periodic extruded mesh. To debug my code / get some understanding, I am typically interested in plotting what's happening on the "axis" $x=y=0$, where $x,y \in D$ are Cartesian coordinates. The simplest way to do this appears to be to use .at() on whatever function and to plot via matplotlib (this may be where I'm going wrong: if there are better ways let me know!).

My issue is that .at() appears to fail for periodic extruded meshes. A minimal working example for me is

from firedrake import *
import numpy as np
import matplotlib.pyplot as plt

# create a flat solid torus
n_tor = 10
refinement_level = 2
periodic = True # .at() works if True, fails if False
disk = UnitDiskMesh(refinement_level = refinement_level) 
mesh = ExtrudedMesh(disk, n_tor, layer_height=2*np.pi/(n_tor), extrusion_type='uniform',periodic=periodic)

# Interpolate a function of the toroidal direction
x,y,z = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, 'Lagrange', 1)
f = Function(V).interpolate(cos(z)) # A periodic function on the torus

# Plot
N = 40*n_tor+1
z_vals = np.linspace(0,2*np.pi,N) + 0.0001 # Bottom contour is sensitive to which side of elements we sample on
pts = np.zeros((N,3))
pts[:,2] = z_vals
f_vals = f.at(pts)

plt.scatter(z_vals, f_vals, label=".at()", color="black")
plt.plot(z_vals, np.cos(z_vals), label="expected")
plt.legend()

This outputs the following plot:

It looks to me like one edge of the elements (in fact, the bottom edge) is evaluating correctly, while the top edges is not. Changing periodic to False completely fixes the problem. I've tried a similar test for SolidTorusMesh to similar results.

Am I going wrong somewhere? Is .at() not available/broken for meshes of this type?

Thanks!

[connorjward]

Does it work if you use a VertexOnlyMesh instead? https://www.firedrakeproject.org/point-evaluation.html#primary-api-interpolation-onto-a-vertex-only-mesh

[maxeruth]

Thanks for the response!

VertexOnlyMesh seems to help. However, there is still a slight issue that I'm running into, which is that VertexOnlyMesh seems to be doing some sort of point reordering when the mesh aligns with cell boundaries. I updated the script to be the following:

from firedrake import *
from firedrake.__future__ import interpolate
import numpy as np
import matplotlib.pyplot as plt

# create a flat solid torus
n_tor = 10
refinement_level = 2
periodic = False # .at() works if True, fails if False
disk = UnitDiskMesh(refinement_level = refinement_level) 
mesh = ExtrudedMesh(disk, n_tor, layer_height=2*np.pi/(n_tor), extrusion_type='uniform',periodic=periodic)

# Interpolate a function of the toroidal direction
x,y,z = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, 'Lagrange', 2)
f = Function(V).interpolate(cos(z)) # A periodic function on the torus

# Interpolate that again to a VertexOnlyMesh
axis_offset = 0.0 # or 0.1
N = 10*n_tor+1
z_vals = np.linspace(0,2*np.pi,N)
pts = np.ones((N,3)) * axis_offset
pts[:,2] = z_vals

vom = VertexOnlyMesh(mesh, pts)
P0DG = FunctionSpace(vom, "DG", 0)
f_vals = assemble(interpolate(f, P0DG))

plt.scatter(z_vals, f_vals.dat.data, label=".at()", color="black")
plt.plot(z_vals, np.cos(z_vals), label="expected")
plt.legend()

The parameter axis_offset is chosen to be 0.0 or 0.1, where 0.0 aligns with the center of the disk. When I use 0.0 I have the following picture:

This doesn't look quite right. However, when I use 0.1, I get the correct picture:

This looks good!

It looks like the problem with the first picture is just that the order of the points got scrambled somehow. I don't know how to recover the correct points however: as far as I can tell this isn't in the linked documentation. I know there is some documentation relating to parallel input ordering, but I am struggling to determine whether that's the problem I am facing (I am running in serial).

[connorjward]

That's very interesting. I don't really know what is causing that. Perhaps the algorithm matching points to cells is faulty. What happens if you refine the mesh?

@leo-collins this might be of interest to you.

[leo-collins]

My guess is also that the cell location algorithm is failing here. With axis_offset=0 your points are exactly on the central vertex of the disk, and for other values of axis_offset your points are on an interior facet. Interesting that it looks to be working with axis_offset=0.1 however; with other values like 0.2, 0.3, 0.4 it looks to be failing again.

If you do something like pts[:,0] *= 0.1 and pts[:,1] *= 0.2 so that your points aren't on an interior facet or vertex it seems to work.

[dham]

The scrambled order is expected. There is no expectation that the order of the vertices in the vertex only mesh is the same as that provided by the user. If you want the values you have created back in the input ordering then you can interpolate them back to a space defined on the input ordering VOM as described in the manual: https://www.firedrakeproject.org/point-evaluation.html#using-the-input-ordering-property

Alternatively, you can work in the vertex only mesh by defining a vector function space and interpolating the vertex only mesh spatial coordinate onto it, hence accessing the locations of the points in VOM order.

[maxeruth]

Hi all, just wanted to say that I got the second option working! Thanks for the help :)