Exact signed distance function of a rhombic dodecahedron

Question

I'm playing around with a ray marcher engine with voxels shaped as rhombic dodecahedra since they can tessallate the 3D space perfectly like hexagons do for the 2D space. The thing is that I only managed to come up with an approximation of the signed distance function for that solid. It's okay to render the voxels but it gives really sharp edges. I would like to make them look better with some rounded edges using an exact SDF.

Does anyone know the exact SDF of a rhombic dodecaheron ?

Here's my (not exact) SDF :

float Sdf(Voxel voxel, vec3 position) {
    position = position - voxel.Position;
    // Exploit the solid's symmetries
    position = vec3(abs(position.x), sign(position.z) * position.y, abs(position.z));

    // distance to each face
    float a = dot(vec3(1.  ,   0.           , 0.          ), position) - voxelSize;
    float b = dot(vec3(0.5 ,   0.           , sqrt3 / 2.  ), position) - voxelSize;
    float c = dot(vec3(0.  , - sqrt2 / sqrt3, 1. / sqrt3  ), position) - voxelSize;
    float d = dot(vec3(0.5 ,   sqrt2 / sqrt3, sqrt3 / 6.  ), position) - voxelSize;
    float e = dot(vec3(0.5 , - sqrt2 / sqrt3, - sqrt3 / 6.), position) - voxelSize;

    return max(a, max(b, max(c, max(d, e))));
}

And this is how it looks like:

EDIT :

So I had some new ideas on how to approach this problem. In this new technique, I try to find the different areas that would have different distance formulas (depending on whether the closest point to the solid is on a face, an edge or a vertex) and compute their distance separately. So far, I only have a correct separation and distance for the areas where the closest point is on a face or on an edge. The code looks like this :

float Sdf(Voxel voxel, vec3 position){
    position = position - voxel.Position;

    vec3 normals[6];
    float dists[6];
    float signs[6];

    // The rhombic dodecahedron has 6 pairs of opposite faces.
    // Assign a normal to each pair of faces.
    normals[0] = vec3(1.  ,   0.           , 0.          );
    normals[1] = vec3(0.5 ,   0.           , sqrt3 / 2.  );
    normals[2] = vec3(0.5 ,   0.           , - sqrt3 / 2.);
    normals[3] = vec3(0.  , - sqrt2 / sqrt3, 1. / sqrt3  );
    normals[4] = vec3(0.5 ,   sqrt2 / sqrt3, sqrt3 / 6.  );
    normals[5] = vec3(0.5 , - sqrt2 / sqrt3, - sqrt3 / 6.);

    // Compute the distance to each face (the sign tells which face of the pair of faces is closest)
    for (int i = 0; i < 6; i++){
        dists[i] = dot(normals[i], position);
        signs[i] = sign(dists[i]);
        dists[i] = max(0., abs(dists[i]) - voxelSize);
    }

    bool sorted = false;
    while(!sorted){
        sorted = true;
        for (int i = 0; i < 5; i++){
            if (dists[i] < dists[i+1]){
                vec3 n = normals[i];
                float d = dists[i];
                float s = signs[i];
                normals[i] = normals[i+1];
                dists[i] = dists[i+1];
                signs[i] = signs[i+1];
                normals[i+1] = n;
                dists[i+1] = d;
                signs[i+1] = s;
                sorted = false;
            }
        }
    }
    
    if (dists[1] < dists[0] * 0.5){ //0.5 comes from sin(pi/6)
        // case where the closest point is on a face
        return dists[0];
    }
    else if (dists[2] == 0.){
        // case where the closest point is on an edge
        float a = 0.5 * dists[0];
        float b = 0.5 * (dists[1] - a);
        return length((dists[0] - b) * signs[0] * normals[0]
                    + (dists[1] - a) * signs[1] * normals[1]);
    }
    else {
        //dunno
    }
}

Here's a quick sketch to show how my approach

The remaining cases get crazier. I'll update the post if any breakthroughs are make.