How to create a Tetrix/Sierpinski Tetrahedron fractal radiating from 0,0,0 ? Python or nodes
Updated answer with animation technique (sort'a cheating)

So to generate an animation of the tetrix levels increasing I kind'a cheated.
Since my previous answer already used a dictionary to store all the levels, I simply wrote that data structure into a json file, then used an app handler function to read it prior to every frame change and load the current level matching the frame number.



So this has solution has two steps:

Generate data file

import numpy as np
from time import time
import bpy, bmesh, json
from mathutils import Vector

def generate_quad( original, quads, level, max_level ):
    ''' Recursive function generating 4 tetrahedra off of the the input tetrahedron '''

    # The 1st new tetrahedron is simply half the scale of the original starting from the original tet's bottom left corner
    origin = original.min(axis=0) / 2    
    half  = original / 2
    half += origin

    # Calculate the dimensions of the new tetrahedra
    newdims = Vector( half.max(axis=0) - half.min(axis=0) )

    # The 2nd moves by the new tet's x dimension to the right
    right = half.copy()
    right[:,0] += newdims.x

    # The 3rd moves by half the x dim and all the y dim
    forward = half.copy()
    forward[:,0] += newdims.x / 2
    forward[:,1] += newdims.y

    # The 4th moves by half the x, a third of the y and all the z
    up = half.copy()
    up[:,0] += newdims.x / 2
    up[:,1] += newdims.y / 3
    up[:,2] += newdims.z

    quad = [ half, right, forward, up ]
    quads[ level ].extend( quad )

    if level &lt; max_level:
        for q in quad:
            quads = generate_quad( q, quads, level + 1, max_level )

    return quads

## Main Code
start = time()
n     = 5  # Number of steps of the fractal
scale = 10 # Scale of the original tetrahedron

# Generate regular tetrahedron at the desired scale
orig = np.array([
    [0,   0, 0],
    [1,   0, 0],
    [0.5, 3**0.5/2, 0],
    [0.5, 1/3*3**0.5/2, ((3**0.5/2)**2 - (1/3*3**0.5/2)**2)**0.5]
]) * scale

# Initialize data structure for all tetrahedra in all steps
quads = { i : [] for i in range( n+1 ) }

# Generate tetrix via recursive function
quads = generate_quad( orig, quads, 0, n )

# From numpy to simple lists to prepare to save as json
quads = { i : np.array(quads[i]).tolist() for i in quads }

quad_cache_file_path = r'C:\Users\username\Documents\tetrix.json'
with open( quad_cache_file_path, 'w' ) as fh:
    quads = json.dump( quads, fh )

print( f'Generated {n} tetrix levels in {round(time() - start,1)} seconds' )


Add frame_change_pre app handler function to link between the frame number and the tetrix's level

import itertools
import numpy as np
import bpy, bmesh, json
from mathutils import Vector

def update_tetrix( scene ):
    C = bpy.context
    D = bpy.data
    S = C.scene

    regular_tet_faces = np.array([[3, 0, 1], [1, 0, 2], [0, 3, 2], [3, 1, 2]])

    quad_cache_file_path = r'C:\Users\username\Documents\tetrix.json'
    with open( quad_cache_file_path ) as fh:
        quads = json.load( fh )

    # Select level to display at current frame (match 
    n = S.frame_current - 1
    n = n if str(n) in quads else max( quads.keys() )

    # Generate mesh data
    verts = []
    faces = regular_tet_faces.tolist()
    for q in quads[str(n)]:
        # Find last vert index and add new verts
        last_vert_index = max([ len( verts ) - 1, 0 ])
        verts.extend( q )

        # Use base face arrangement of original regular tetrahedron and offset indices
        faces_indices = regular_tet_faces + last_vert_index + 1
        faces.extend( faces_indices.tolist() )

    # Remove duplicates
    faces = sorted([ sorted(f) for f in faces])
    faces = list(f for f,_ in itertools.groupby(faces))

    if 'tetrix' in D.objects:
        print( f"Updating mesh with level {n}" )
        o = D.objects['tetrix']
        m = o.data

        bm = bmesh.new()    

        bmverts = [ bm.verts.new(co) for co in verts ]
        bmfaces = [ bm.faces.new([ bmverts[i] for i in f]) for f in faces ]

        bm.to_mesh( m )    
        m.update()       

    else:
        # Generate mesh data and initialize vertex and face arrays
        m = D.meshes.new('tetrix')

        # Generate mesh object and add to scene
        m.from_pydata( verts, [], faces )
        print({ k : len( quads[k] ) for k in quads })
        o = D.objects.new('tetrix',m)
        S.collection.objects.link( o )

bpy.app.handlers.frame_change_pre.append(update_tetrix)


Caveats


This solution doesn't enable animating a transition between levels, simply a harsh switch. You also don't get a lot of control over how long the pauses are between switches (though this could be done in post in any video editing software by multiplying or delaying frames).
Currently the entire data file is read prior to every frame change. This can be pretty slow if you've generate a big file with many levels (above 7), depending on your computer's specs. It could probably be improved by loading the file once to memory as a global scene object then reading directly from that prior to every frame change, rather than reading the file from disk every time.


Original Answer



This script can generate a n-level tetrix. It will work only in Blender 2.8 and utilizes some Numpy tricks for more efficient numerical manipulation.

import bpy, bmesh
import numpy as np
from mathutils import Vector

def generate_quad( original, quads, level, max_level ):
    ''' Recursive function generating 4 tetrahedra off of the the input tetrahedron '''

    # The 1st new tetrahedron is simply half the scale of the original starting from the original tet's bottom left corner
    origin = original.min(axis=0) / 2    
    half  = original / 2
    half += origin

    # Calculate the dimensions of the new tetrahedra
    newdims = Vector( half.max(axis=0) - half.min(axis=0) )

    # The 2nd moves by the new tet's x dimension to the right
    right = half.copy()
    right[:,0] += newdims.x

    # The 3rd moves by half the x dim and all the y dim
    forward = half.copy()
    forward[:,0] += newdims.x / 2
    forward[:,1] += newdims.y

    # The 4th moves by half the x, a third of the y and all the z
    up = half.copy()
    up[:,0] += newdims.x / 2
    up[:,1] += newdims.y / 3
    up[:,2] += newdims.z

    quad = [ half, right, forward, up ]
    quads[ level ].extend( quad )

    if level &lt; max_level:
        for q in quad:
            quads = generate_quad( q, quads, level + 1, max_level )

    return quads

## Main Code
C = bpy.context
D = bpy.data
S = C.scene

n     = 2  # Number of steps of the fractal
scale = 10 # Scale of the original tetrahedron

# Generate regular tetrahedron at the desired scale
orig = np.array([
    [0,   0, 0],
    [1,   0, 0],
    [0.5, 3**0.5/2, 0],
    [0.5, 1/3*3**0.5/2, ((3**0.5/2)**2 - (1/3*3**0.5/2)**2)**0.5]
]) * scale

regular_tet_faces = np.array([[3, 0, 1], [1, 0, 2], [0, 3, 2], [3, 1, 2]])

# Initialize data structure for all tetrahedra in all steps
quads = { i : [] for i in range( n+1 ) }

# Generate tetrix via recursive function
quads = generate_quad( orig, quads, 0, n )

# Generate mesh data and initialize vertex and face arrays
m     = D.meshes.new('tetrix')
verts = []
faces = regular_tet_faces.tolist()

# Generate mesh data
for q in quads[n]:
    # Find last vert index and add new verts
    last_vert_index = max([ len( verts ) - 1, 0 ])
    verts.extend( q.tolist() )

    # Use base face arrangement of original regular tetrahedron and offset indices
    faces_indices = regular_tet_faces + last_vert_index + 1
    faces.extend( faces_indices.tolist() )

# Generate mesh object and add to scene
m.from_pydata( verts, [], faces )
print({ k : len( quads[k] ) for k in quads })
o = D.objects.new('base',m)
S.collection.objects.link( o )


Here's a higher level tetrix from all orthographic angles and a free one (quad view):