How to create a Tetrix/Sierpinski Tetrahedron fractal radiating from 0,0,0 ? Python or nodes

Question

Animation nodes method true to description has been provided by @lemon, Vote his answer UP UP UP!

INFERIOR

INFERIOR

INFERIOR / different concept

Mesh. I want to create it not by dissection of edges, but instead by splitting 1 into 4, or should I say, adding 3 more to the original location before moving all 4 away from Local 0,0,0 , and then doing the same with the 4 replicas, as each of them is splitting the same way (the subsequent movements of replicas becoming more and more complex/derived), so that (if animated) it really looks like the one is splitting into 4, and each of the 4 results are then splitting into 4 (and so on).
If Empties were placed in each start/end position of each object it would, looking at the Empties hierarchically connected, appear like a branching tree, with old branches continuously lengthening.

Though I wouldn't mind knowing the opposite too. Having 1 edges dissected(halved) and created vertexes connected, but always so that connection doesn't happen toward the inside (coordinate 0,0,0), THOUGH THAT WOULDN'T really be the opposite of the previous, the real opposite would look like the 1 is scaling in (previously mentioned) 4 directions (toward the vertexes of the 1), and so on, in each of the 4 newly created. Yeah XD definitely wouldn't mind knowing how to do this.

Method using python script or nodes would be just FINE!

IRL construction

BTW I have a hunch it would be easier to construct it if there was some alignment with the coordinate system, such as edge midpoints aligning on the axes.

Answer 1

Animation nodes version

One tetrahedron is used as model.

Another one is used as duplivert. Its parent's vertices will be constructed by the nodes.

An empty is used to control the iterations. And this is triggered when the empty is scaled.

The nodes setup is composed of 4 parts.

• A main: get the input and produce the resulting vertices into the dupli parent
• A first sub part: calculate the scale factors and iterate to construct each level
• A second sub part: iterates on the model to produce each vertex
• And a last sub part: duplicate each vertex of the model and place it

The scales are calculated by power of two: 0, 1, 2, 4, etc. For each level, depending on the scale of the empty, we have a scale that is proportional to the intervals [0,1], [1,2], [2, 4], etc. The final unit is the base tetrahedron obtained from regular solid addon.

Main

This main part calculates the amount of iteration and the relative scale. For instance, if the empty is scale by 3.20, we have 3 iterations and a position of 0.2 for the scale.

First sub part

From that the iteration is done from a base vertex input which is reassigned to itself after each iteration.

The calculation is to determinate the scale level ([0,1], [1,2], [2, 4], etc.) and adjust it from the scale position (the 0.2 above).

Second and third sub parts

Simply iterates on the model for each cumulated vertex.

And in the last sub part, place the vertex scaled and offset from the model.

When the empty is scaled up to ten, that will generate more than 5 millions vertices. So... be careful.

Python version

Uses the same principles as the AN version above.

I've added the possibility to use other shapes.

So we can have:

Or:

The principle is to use:

• An empty to scale the figure (scale the empty to scale it)
• A model object which gives the base vertices
• A dupli object which gives the shape of the inner elements
• A target object (dupli parent of the dupli) for which the code calculates the vertices

The script itself is in the following parts:

• Setup a hook function on bpy.app.handlers.scene_update_pre so that it can be animated
• Get the key objects (hard coded) to apply the generative code
• The generation itself
• A function to clean the target mesh
• A function to calculate the wanted scales depending of the depth of the figure

Scene handler

Remove previous ones and setup/update the hook

def setup_handler():
    for h in [x for x in bpy.app.handlers.scene_update_pre]:
        if h == update_sierpinsky:
            bpy.app.handlers.scene_update_pre.remove( update_sierpinsky )

    bpy.app.handlers.scene_update_pre.append( update_sierpinsky )

Get key objects

Here if you change (in the given .blend below) model1/dupli1 by model2/dupli2 you'll have cubes instead of Sierpinski.

def update_sierpinsky( scene ):
    scale_object = scene.objects['scale']
    model_object = scene.objects['model1']
    dupli_object = scene.objects['dupli1']
    target_object = scene.objects['target']

    sierpinsky( scale_object, model_object, dupli_object, target_object )

Scales calculations

Calculate 0, 1, 2, 4, etc. plus a value between 0 and 1 to make the transitions. More precisely, this is 0, 1, sqrt of the model vertex count as the scales are not the same depending on the figure (but not sure of this formula...).

def scales( scale_object, scale_factor ):
    scale = scale_object.scale.z - 1
    scale = max( scale, 0 )
    scale_floor = int( scale )
    scale_delta = scale - scale_floor
    iterations = scale_floor + 1

    def calc( sd, i, sf ):
        i2 = sf ** i
        prev = i - 1
        prev2 = sf ** prev if i != 0 else 0
        return prev2 + sd * (i2 - prev2)

    return [calc( scale_delta, i, scale_factor ) for i in reversed( range( iterations ) )]

Vertices generation

I've choose an iterative one (for a variant).

It loops overs the scales and generates the vertices locations level by level.

def sierpinsky( scale_object, model_object, dupli_object, target_object ):

    model_vertices = [Vector( v.co ) for v in model_object.data.vertices]
    result_vertices = [Vector()]
    scale_factor = len(model_vertices) ** 0.5

    for scale in scales( scale_object, scale_factor ):
        scaled_vertices = [scale * v for v in model_vertices]
        result_vertices = [r + s for r in result_vertices for s in scaled_vertices]

    clear_mesh( target_object )
    target_object.data.from_pydata( result_vertices, [], [] )
    target_object.data.update()

Target mesh clean up

The vertices of the target mesh are removed before each calculation.

There may be a better way to do that but don't know how.

def clear_mesh( obj ):
    bm = bmesh.new()
    bm.from_mesh( obj.data )
    bm.verts.ensure_lookup_table()
    for v in [x for x in bm.verts]: # use a copy before removing
        bm.verts.remove(v)
    bm.to_mesh( obj.data )
    obj.data.update()

The blend file for the Python version:

Answer 2

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 < 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

1. 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).

2. 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 < 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):

Answer 3

Bmesh version.

Pyramid

Scaling a tetra particle on a single vert point cloud, such that each time it doubles in size, it uses the next level pyramid as a particle

Make one tetrahedron from coordinates. I have skinned it lazilly with convex hull

Now for each level, its pretty much from any vertex corner of the first tetrahedron, make 3 copies and move to the other end of each linked edge.

As the level increases by one the number of unit tetra edges that make up the offset doubles.

So pretty much make a unit mesh, for each level copy and distribute mesh 3 times, write to mesh.

Pointcloud

As well as creating a skinned unit form tetrahedrons, can simply seed the pyramid with a single vertex mesh, at the tetra origin. Dupliverts or particle systems can then instance our units.

Four units, a level 0, 1, 2 and 3 pyramid. Three point clouds at level 0, 4, and 8

Script to make the 4 units and 3 point clouds.

import bpy
import bmesh
from mathutils import Vector
context = bpy.context
scene = context.scene

levels = 12
sublevels = 4


def new_coll(name, parent=None, obs=()):
    c = bpy.data.collections.new(name)
    if parent:
        parent.children.link(c)
    for o in obs:
        c.objects.link(o)
        o.hide_set(True)
    return c


coords = [Vector(c) for c in (
    (0, 0, 0),
    (1, 0, 0),
    (.5, 3**.5 / 2, 0),
    (.5, 1 / 3 * 3**.5 / 2, ((3**.5 / 2)**2 - (1 / 3 * 3**.5 / 2)**2)**.5))]

o = sum(coords, Vector()) / 4


def pyramid(me, levels, origin=o, scale=(1, 1, 1)):
    bm = bmesh.new()
    bm.from_mesh(me)
    for l in range(1, levels + 1):
        lvrts = len(me.vertices)
        for vec in coords[1:]:
            bm.from_mesh(me)
            bmesh.ops.translate(bm,
                                vec=2 ** (l - 1) * vec,
                                verts=bm.verts[-lvrts:],
                                )
        bm.to_mesh(me)
        # me.update()

    # to change origin to centroid.
    x = sum([v.co for v in bm.verts], Vector()) / len(bm.verts)

    bmesh.ops.translate(bm,
                        vec=-x,
                        verts=bm.verts,
                        )

    bmesh.ops.scale(bm,
                    vec=scale,
                    verts=bm.verts)

    bm.to_mesh(me)
    bm.free()
    return me


bm = bmesh.new()
# make a tetra
bmesh.ops.convex_hull(bm,
                      input=[bm.verts.new(c - o) for c in coords]
                      )
me = bpy.data.meshes.new("Tetra")
bm.to_mesh(me)
ob = bpy.data.objects.new("0Tetra", me)


trap_coll = new_coll("Trapinski", scene.collection)
units = [ob]
for level in range(1, sublevels):
    units.append(bpy.data.objects.new(f"{level}Pyramid",
                                      pyramid(me.copy(), level, origin=o / 4, scale=(2 ** -(level), ) * 3)))

unit_coll = new_coll("Units", trap_coll, units)

bm.clear()
frames = []
for level in range(0, levels, sublevels):  # levels):
    mef = bpy.data.meshes.new(f"Frame{level}")
    bm.verts.new(o)
    bm.to_mesh(mef)
    pyramid(mef, level)
    ob = bpy.data.objects.new(f"Frame{level}", mef)
    ob.scale = (2 ** -level,) * 3
    frames.append(ob)
    ob.instance_type = 'VERTS'


frame_coll = new_coll("Frames", trap_coll, frames)

Animating

Using a method similar to @TLousky with a frame change handler. continuing from script above

Swaps pointcloud every 40 frames and unit parented to that pointcloud ever 10 frames.

def show(units, frames):
    def handler(scene):
        f = scene.frame_current
        uidx = (f % 40) // 10
        idx = (f // 40) % len(frames.objects)
        #uidx = 0
        #uidx = scene.get("prop", 0)
        unit = units.objects[uidx]
        for frame in frames.objects:
            frame.hide_set(True)
        frame = frames.objects[idx]
        frame.hide_set(False)

        for u in units.objects:
            u.scale = (1, 1, 1)
            u.parent = None
        unit.parent = frame
        frame.hide_set(False)
    return handler


bpy.app.handlers.frame_change_pre.clear()
bpy.app.handlers.frame_change_pre.append(show(unit_coll, frame_coll))

Particle System

Another way of instancing is via a particle system.

Here I've added a simple particle system to single vertex frame 0 and shown result using each of 4 units as particle object

Can stack a number of particle systems per pointcloud objects, staggering a start and end frame and life time.

Using particles to grow the tree.

Growing the pyramid like a tree

AFAIC this is a tree structure. It doesn't radiate out from its centre, rather it grows from the selected base.

Stepping through the first few frames. The top vertex is vertex group "level0" and particle system "level0" of the point cloud, the tetra particle grows (particle scale 0 to 1) in 5 frames. (given a long life time or they start vanishing (cool effect actually), whereapon the next level is used, and so on

Here I've made one tetra, and one point cloud. Each level of the cloud is assigned a vertex group, and for each vertex group used as density for a staggered frame start particle system. The particle is keyframed to scale from 0 to 1 over tpf frames. For a level levels point cloud.

import bpy
import bmesh
from mathutils import Vector
context = bpy.context
scene = context.scene


def new_coll(name, parent=None, obs=()):
    c = bpy.data.collections.new(name)
    if parent:
        parent.children.link(c)
    for o in obs:
        c.objects.link(o)
        #o.hide_set(True)
    return c


coords = [Vector(c) for c in (
    (0, 0, 0),
    (1, 0, 0),
    (.5, 3**.5 / 2, 0),
    (.5, 1 / 3 * 3**.5 / 2, ((3**.5 / 2)**2 - (1 / 3 * 3**.5 / 2)**2)**.5))]

o = coords[3]


def pyramid(me, levels, origin=None, scale=(1, 1, 1)):
    bm = bmesh.new()
    bm.from_mesh(me)
    for l in range(1, levels + 1):
        lvrts = len(me.vertices)
        for vec in coords[1:]:
            bm.from_mesh(me)
            bmesh.ops.translate(bm,
                                vec=2 ** (l - 1) * vec,
                                verts=bm.verts[-lvrts:],
                                )
        bm.to_mesh(me)
        # me.update()

    # to change origin to centroid.
    if origin is not None:
        x = Vector(o)
    else:
        x = sum([v.co for v in bm.verts], Vector()) / len(bm.verts)

    bmesh.ops.translate(bm,
                        vec=-x,
                        verts=bm.verts,
                        )

    bmesh.ops.scale(bm,
                    vec=scale,
                    verts=bm.verts)

    bm.to_mesh(me)
    bm.free()
    return me


bm = bmesh.new()
# make a tetra
bmesh.ops.convex_hull(bm,
                      input=[bm.verts.new(c - o) for c in coords]
                      )
me = bpy.data.meshes.new("Tetra")
bm.to_mesh(me)
ob = bpy.data.objects.new("0Tetra", me)
# parameters ##########################
fpt = 5
level = 6
#######################################
trap_coll = new_coll("Trapinski", scene.collection)
units = [ob]
unit_coll = new_coll("Units", trap_coll, units)
bm.clear()
frames = []
ob.hide_set(True)
mef = bpy.data.meshes.new(f"Frame{level}")
bm.verts.new()
bm.to_mesh(mef)
pyramid(mef, level)
f = bpy.data.objects.new(f"Frame{level}", mef)

frames.append(f)

bm.clear()
bm.from_mesh(mef)
verts = sorted(bm.verts, key=lambda v: v.co.z)

lev = 0
while verts:    
    v = verts.pop()
    print(v.co)
    vlevel = [v.index]
    while verts:        
        if abs(v.co.z - verts[-1].co.z) < 1e-6:
            vlevel.append(verts.pop().index)
        else:
            break
    vg = f.vertex_groups.new(name = f"level{0}")
    vg.add(vlevel, 1, 'REPLACE')
    ps = f.modifiers.new(f"level{lev}", 'PARTICLE_SYSTEM')
    psys = ps.particle_system
    psys.settings.count = len(vlevel)
    psys.settings.frame_start = lev * fpt + 1
    psys.settings.frame_end = lev * fpt + 1
    psys.settings.lifetime = 1000
    psys.settings.render_type = 'OBJECT'
    psys.settings.instance_object = ob
    psys.settings.emit_from = 'VERT'
    psys.settings.use_emit_random = False
    psys.settings.physics_type = 'NO'
    psys.settings.particle_size = 0
    psys.settings.keyframe_insert("particle_size", frame=lev * fpt + 1)
    psys.settings.particle_size = 1
    psys.settings.keyframe_insert("particle_size", frame=lev * fpt + fpt + 1)
    psys.settings.use_global_instance = True
    psys.settings.use_rotation_instance = True
    psys.settings.use_scale_instance = True

    psys.settings.rotation_mode = 'NONE'
    psys.vertex_group_density = vg.name

    lev += 1    


frame_coll = new_coll("Frames", trap_coll, frames)

Note could use a certain level staggered to be the next level, ie use the animation created again offset to the bottom three corners.

Collection Instance

Lastly, can make a collection instance from the duplivert frame / unit pair.

The collection instance of frame4 dupliverting simple tetra

This instance can be used as a particle.

Scaling

The code above is set up to have the units unit size, and the frames scaled such that the whole is unit size. The pyramid will double in size for each level, and gets large very quickly. It is not splitting edges.

Adding a parent empty to double unit in size on every level change. Looking from above in viewport, till it gets so large we are inside it looking thru the bottom gap, at which time the viewport is rolled around x

Summing Up.

There are a number of techniques shown here to make animations going from base tetra to extremely high level pyramids, only using one tetra, or a few simple levels of faced meshes.

Answer 4

I consider TLousky's answer superior, especially for pedagogic reasons (expressiveness and documentation of the code).

A few minutes late but without numpy (also 2.8 only):

import bpy
from itertools import combinations
from mathutils import *

scale = .333 #scale of one tetrahedon
tetrons = 0  #index of current tetron
steps = 4    # recursion steps

vec =(
(0,0,0),
(1,0,0),
(.5,3**.5/2,0),
(.5,1/3*3**.5/2, ((3**.5/2)**2 - (1/3*3**.5/2)**2)**.5))
verts = Matrix(vec)
vertices = []
faces = []

def Tetron(vec=Vector([0,0,0]), ):
    #Adds tetrahedons to vertices and faces
    global tetrons
    base_vecs = verts*scale
    for base_vec in base_vecs:
        base_vec += vec*scale
        vertices.append(base_vec)
    faces_ = list(combinations(list(range(tetrons*4, tetrons*4+4)), r = 3))
    faces.extend(faces_)
    print("ind {}".format(tetrons))
    tetrons+=1

def SiPyramid(depth, vec=Vector([0,0,0])):
    #Recursive function to add Tetrahedons
    if depth < 0:
        return 0 
    if depth==0:
        Tetron(vec)
    else:
        for base_vec in verts:
            print(depth-1)
            print((2**(depth-1)*base_vec + vec ))
            SiPyramid(depth-1, 2**(depth-1)*base_vec + vec   )

SiPyramid(steps)

#Link data to blender objects
mesh_data = bpy.data.meshes.new("sierpinsky")
mesh_data.from_pydata(vertices, [], faces)
mesh_data.update()
obj = bpy.data.objects.new("Sierpinsky", mesh_data)
bpy.context.collection.objects.link(obj)

At 10 steps (might take a few seconds to generate):

Answer 5

I've produced a Sierpinski Pyramid using Array modifiers. Two are needed per iteration.

First, create a tetrahedron. This can be done by creating a cube, deleting two opposite vertices from the top face and the other two opposite vertices from the bottom face, then reconnecting the remaining vertices. Then rotate it and align it so that one face is flat and the top vertex is the object's origin.

Next create an empty at the origin, and rotate it 120 degrees. Call this RotEmpty.

Finally create an empty at one of the bottom vertices. Call this LocEmpty1.

Apply two Array modifiers to the object. The first should use LocEmpty1 for its offset and have a count of 2. The second should use RotEmpty and have a count of 3. Check the boxes to merge vertices - this will help reduce the geometry as it's not optimal...

To add another tier, place a new Empty at the bottom of the pyramid, duplicate both modifiers on the pyramid and set the Loc one to the new empty. Repeat as needed until you have the number of tiers you need.

The advantage to generating it this way is that you can animate it. I did this by parenting everything to the pyramid, then animating each Empty going from the origin to its "full" location in sequence, and moving the whole pyramid up to keep its center at the origin. Set the camera to zoom out over the whole thing, slap a sun light on there, and... ta-da!

Answer 6

This is a non-script approach. You can create this fractal adding a tetrahedron with Extra Objects add-on. Then subdivide it once, remove middle triangles and make inner faces. Repeat the subdivide once and remove and make faces steps a few times. Scale the fractal down and then animate scaling it up to show the structure.