Thanks to this question and answers I could put together this code (The configurations.h5 file comes from the known solutions of the Thompson problem).

If you execute the code blender3 generates the base shape for a d9. However I still miss the final steps, would be possible to add them in python?

Here how to do manually:

Go under Modelling, select the mesh, go in Edit mode, and then you can add the missing faces: press 3af (3 to select faces mode, a for select all, f to add faces).

Now, for all the points in the variable points inverted we need to add the numbers, as example the opposite of the pole made by the code (0, 0,  -1).

Under Modelling, add Text, set up coordinates (0, 0, 1), put Alignment center both in Vertical and Horizontal, fix font and size. Convert to Mesh, and finally in edit mode extrude inside the die.

To dig the text from the die you can finally use the Difference Boolean tool: remove the extruded text from the die.

Now, I need to repeat the text for every value. While doable for a d9, it would be extremely demanding for larger dice, at very least to copy correctly coordinates and angles. So, as mentioned would be possible to put those steps in the python code? How?

import bpy
import bmesh
from mathutils import Vector
import h5py

n = 9 #n is the number of points
rat = .86 # how far along the radius to bisect

u_segments = 128  # UV sphere settings
v_segments = 128

h5 = h5py.File('/full/path/to/configurations.h5','r')
hset = h5[str(n)]
coords = hset.get('coordinates')
ocoords = []
for item in coords:
    ocoords.append(item[0])
    ocoords.append(item[1])
    ocoords.append(item[2])
h5 = None
hset = None
coords = None

it = iter(ocoords)
points = list(map(Vector, zip(it, it, it)))
ocoords = None

context = bpy.context
scene = context.scene
# make one point south pole.
R = points[0].rotation_difference(Vector((0, 0, -1))).to_matrix()
points = [R @ p for p in points]

bm = bmesh.new()
bmesh.ops.create_uvsphere(bm,
        radius=1,
        u_segments=u_segments,
        v_segments=v_segments)

for p in points:
    ret = bmesh.ops.bisect_plane(bm,
            geom=bm.faces[:]+bm.edges[:]+bm.verts[:],
            plane_co= rat * p,
            plane_no=-p,
            clear_outer=False,
            clear_inner=True)

me = bpy.data.meshes.new("dice")
bm.to_mesh(me)
ob = bpy.data.objects.new("dice", me)
scene.collection.objects.link(ob)
context.view_layer.objects.active = ob
ob.select_set(True)
ob.location = scene.cursor.location

The generic die is build with a clever idea: to flatten N equally distributed circles in a sphere. .

Throwing the chopped sphere will land on one of the flattened parts and you should be able to read the result on top.

The shape is a bit weird, but thanks to this question and answers I could put together the code you can see in the bottom (The configurations.h5 file comes from the known solutions of the Thompson problem).

If you execute the code blender3 generates the base shape for a d9. This one:

The shape is a great start, but I still miss the final steps. I'm doing the final steps in GUI, how to do it in Python automatically?

Here how to do in GUI:

First I need to close the holes. Go under Modelling, select the mesh, go in Edit mode, and then you can add the missing faces: press 3af (3 to select faces mode, a for select all, f to add faces).

The mesh become like this:

After covering the hole I need to put the digits. The digits have to appear in the opposite side compared to the flatten circles.

A flatten circle is always in the bottom (coord 0,0,-1), so let's start from that coordinate for it is easier.

Under Modelling, add Text, set up coordinates (0, 0, 1), put Alignment center both in Vertical and Horizontal, fix font and size.

Convert to Mesh, and finally extrude and move inside the die.

To dig the text from the die you can finally use the Difference Boolean tool: remove the extruded text from the die.

Now, I need to repeat the text for every value. While doable for a d9, it would be extremely demanding for a larger dice, at very least to copy correctly coordinates and angles (in the example are 0,0,1 and 0,0,0. But other faces would have other numbers.)

So, would be possible to put those steps in the python code? How?

import bpy
import bmesh
from mathutils import Vector
import h5py

n = 9 #n is the number of faces
rat = .86 # how far along the radius to bisect

u_segments = 128  # UV sphere settings
v_segments = 128

h5 = h5py.File('/full/path/to/configurations.h5','r')
hset = h5[str(n)]
coords = hset.get('coordinates')
ocoords = []
for item in coords:
    ocoords.append(item[0])
    ocoords.append(item[1])
    ocoords.append(item[2])
h5 = None
hset = None
coords = None

it = iter(ocoords)
points = list(map(Vector, zip(it, it, it)))
ocoords = None 


context = bpy.context
scene = context.scene
# make one point south pole.
R = points[0].rotation_difference(Vector((0, 0, -1))).to_matrix()
points = [R @ p for p in points] 


bm = bmesh.new() 

bmesh.ops.create_uvsphere(bm,
        radius=1,
        u_segments=u_segments,
        v_segments=v_segments)

for p in points:
    ret = bmesh.ops.bisect_plane(bm,
            geom=bm.faces[:]+bm.edges[:]+bm.verts[:],
            plane_co= rat * p,
            plane_no=-p,
            clear_outer=False,
            clear_inner=True)

me = bpy.data.meshes.new("dice")
bm.to_mesh(me)
ob = bpy.data.objects.new("dice", me)
scene.collection.objects.link(ob)
context.view_layer.objects.active = ob
ob.select_set(True)
ob.location = scene.cursor.location