FEniCS, pygmsh and boundary conditions

Well, after all, I think it is due to the boolean operation:

sphere_inner = geom.add_ball([sphere_x1, sphere_y1, sphere_z1], sphere_r1)
sphere_outer = geom.add_ball([sphere_x2, sphere_y2, sphere_z2], sphere_r2)
sphere_final = geom.boolean_difference([sphere_outer],[sphere_inner], delete_first=True, delete_other=True)
geom.add_raw_code("Physical Surface(\"1\") = {1,2,3};") # inner
geom.add_raw_code("Physical Surface(\"2\") = {4,5,6};") # outer
geom.add_raw_code("Physical Volume(\"3\") = {" + sphere_final.id + "};") # volume

I subtract the small from the large sphere, which results into a sphere with a spherical ‘hole’ (this is what I expect). Then I assign the surfaces and volume.

However, Gmsh tells me now that the inner sphere surface has the same Physical surface number as the outer one (see image), which should not be … . I will see the next days … .