Thanks for your answer. I tried what you posted. Seems that everything works but I am not sure about something, that is probably trivial. Consider, for example, a simple cone:

Im not sure what you expect.
You are only saving triangles (which meshio only picks those marked as entity groups in gmsh).
If you want to mesh the 3D mesh, I would expect you to save tetrahedral cells, i.e.

I have a hopefully final question, thankfully not very naive. In my legacy fenics code I use for a variety of reasons external functions for my boundaries that have numerical calculations on them. I did something like this

class BCfunction(UserExpression):
def eval(self, value, x):
value[0] = f(x) # some other extrernal function of x with non trivial numerics
def value_shape(self):
return ()
u_D = BCfunction()
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)

However I am not sure what is the better way to do this in fenicsx

You would be more explicit in DOLFINx.
You should:

create a function u_bc = dolfinx.fem.Function(V).

Create a python function

def bc_func(x):
return f(x)

f(x) is here assumed to be vectorized, i.e. it takes in a 2D array of points x = ((x_0, x_1, ... x_N), (y_0, y_1, ..., y_N), (z_0, z_1, ...., z_N)
If your function is not vectorized, you would need to transpose the input

def bc_func(x):
xT = x.T
values = np.zeros(x.shape[1])
for (i,coord) in enumerate(xT)):
values[i] = f(coord)
return values

As a last note, I guess if my function f(x) calculates the boundary conditions only at the actual boundaries of my geometry (that is why I asked how to get them in my original question), i.e. an analytical function in position space x but with parameters from other calculations, is the interpolation step optional?