Could test cylindrical and spherical exports with more complex geometries

[jhdark]

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:

@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)