Could test cylindrical and spherical exports with more complex geometries

We should also consider geometries which also vary in r and z in cylindrical and spherical coordinate systems

here is an example test that could be adapted:
```python
@pytest.mark.parametrize("length", [4, 5, 3])
def test_average_surface_rhombus(length):

# creating a mesh with FEniCS
# Define the number of divisions in the x and y directions
nx, ny = 10, 10
mesh = f.RectangleMesh(f.Point(0, 0), f.Point(length, length), nx, ny)

# Access mesh coordinates
coordinates = mesh.coordinates()

# Rotation matrix for 45 degrees
theta = np.pi / 4 # 45 degrees in radians
rotation_matrix = np.array(
[[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]
)

# Apply rotation to each coordinate
for i, coord in enumerate(coordinates):
# Rotate the point (x, y) using the rotation matrix
coordinates[i] = np.dot(rotation_matrix, np.array(coord))

surface_markers = f.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
surface_markers.set_all(0)

# find all facets along y = x and mark them
surf_id = 2
for facet in f.facets(mesh):
midpoint = facet.midpoint()
if np.isclose(midpoint[0], midpoint[1], atol=1e-10):
surface_markers[facet.index()] = surf_id

ds = f.Measure("ds", domain=mesh, subdomain_data=surface_markers)
my_export = AverageSurface("solute", surf_id)
V = f.FunctionSpace(mesh, "P", 1)

c_fun = lambda x, y: 2 * x + y
expr = f.Expression(
ccode(c_fun(x, y)),
degree=1,
)
my_export.function = f.interpolate(expr, V)
my_export.ds = ds

expected_value = (3 * length) / (2 * (2**0.5))
computed_value = my_export.compute()

assert np.isclose(computed_value, expected_value)
```