Partitioning a Cartesian Product Yields Sierpiński Triangle... ish?

Question

I wanted to find an efficient method to partition a Cartesian product of $n$ sets $S_i$ of varying sizes into maximum size subsets that are defined by all tuples in the partition differing in at least $k$ places. I failed to come up with a non-bruteforce algorithm. But I was pretty sure there should be a better approach. So I computed a few partitions by brute force and visualized them to maybe get a new idea. I found that for the case where $k = 2$ using 3 sets of size 100 each, the result seems to be something like what is apparently called tetrix, a 3-D Sierpiński Triangle.

I don't understand this result and it did not help me in coming up with an algorithm. But since this is a fractal I am now even more sure there needs to be a non-bruteforce computation to generate this fractal and maybe also a generalization for all possible setups of sets $S_i$ and $k$.

Can anyone explain the result to me and maybe knows how to non-brutoforce generate these partitions?

This is the code:

import seaborn as sns, numpy as np, pandas as pd, random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from tqdm import tqdm

L = [list(range(100)), list(range(100)), list(range(100))]

from itertools import product
def n_distinct_product(iterables, n_distinct):
    elements_n_distinct = set()
    pr = list(product(*iterables))
    for s1 in tqdm(pr):
        if all(sum([w1 != w2 for w1, w2 in zip(s1, s2)]) >= n_distinct
               for s2 in elements_n_distinct if s1 != s2):
            elements_n_distinct.add(tuple(s1))
    return elements_n_distinct


k = 2
data = n_distinct_product(L, k)
data = list(data)

data = np.vstack(data)

df = pd.DataFrame(data=data, columns = ['x', 'y', 'z'])


sns.set_style("whitegrid", {'axes.grid' : False})

fig = plt.figure(figsize=(6,6))

ax = Axes3D(fig)


g = ax.scatter(df.x, df.y, df.z, c=df.z, marker='.', depthshade=True, cmap='Paired')
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

# produce a legend with the unique colors from the scatter
legend = ax.legend(*g.legend_elements(), loc="lower center", title="X Values", borderaxespad=-10, ncol=4)
ax.add_artist(legend)

plt.show()

sns.pairplot(data=df)

sns.scatterplot('x', 'y', data=df, hue='z')

Description of algorithm:

Define a set $M = \{\}$. Iterate over Cartesian product in lexicographic order, starting with $(0,0,0)$. For each tuple: If it differs from all tuples in $M$ in at least $k$ places, add it to $M$.