How to edit boundary conditions from python

Hello everybody,
that would be great if someone could give me a hint on the following. I need to assign boundary conditions to the faces of my mesh, exactly what is done by Mesh->Edit Boundary Condition... in NGSolve GUI, but from a script, because I will have to do it many times and there are many faces.

I am trying to do it with a python script, a simplified version of which looks like this:

from netgen.meshing import *
from netgen.csg import *
from ngsolve import *

geometry = CSGeometry()

bigball = Sphere(Pnt(0.0, 0.0, 0.0), 100.0);
smallball = Sphere(Pnt(0.0, 0.0, 0.0), 20.0);

vac = bigball - smallball

geometry.Add(vac)

ngmesh = geometry.GenerateMesh()

ngmesh.Add(FaceDescriptor(surfnr=1,domin=1,bc=2))
ngmesh.Add(FaceDescriptor(surfnr=2,domin=1,bc=1))

ngmesh.Export('.', 'Elmer Format')

mesh = Mesh(ngmesh)

Draw(mesh)

This script can be executed, but it will not assign boundary 2 to face 1 and boundary 1 to face 2, as I would like it to. Instead, it will assign:
face -- boundary
1 -- 1
2 -- 2
3 -- 2
4 -- 1
So it simply adds non-existing faces 3 and 4 instead of changing the assignments for faces 1 and 2. Does anybody know how to make netgen change boundary numbers assigned to faces instead of adding new faces?

Or any other way to change these assignments automatically? I thought about replacing the boundary numbers in the mesh.boundary file afterwards, but netgen only assigns boundary numbers from 1 to 22 and I have more faces, so I can't identify faces by boundary numbers.

You can give the bc numbers to the sphere object:
Best
Christopher

Thank you Christopher!
It has worked with spheres, however for more complicated cases it seems to not work properly. In the following example only two out of six faces of the second brick get boundary condition 3, the rest are 2. What could be the problem here?

I didn't manage to find out how this boundary condition assignment works if constructive solid geometry (CSG) and boolean operations are used. I have just written a python script that rewrites boundary condition assignments to every mesh element in mesh.boundary file according to the nodes' coordinates from the mesh.nodes file. Thus I have achieved what I wanted.