The following codes show examples to create polygon mesh geometries with branching algorithm. The first code below shows the basic branching algorithm on 2D as a start.

add_library('igeo')

def setup() : 
    size(480, 360, IG.GL)
    IG.duration(120)
    MyAgent(IG.v(0, 0, 0), IG.v(1, 0, 0)).clr(0)

class MyAgent(IAgent) : 
    def __init__(self, pos, dir) : 
        self.pt1 = pos
        self.pt2 = pos.cp(dir)
        self.isColliding = False
    
    def interact(self, agents) : 
        for agent in agents : 
            if self.isColliding :
                return
            if isinstance(agent, MyAgent) : 
                if agent is not self : 
                    dist = agent.pt2.dist(self.pt2) #distance of end points
                    tolerance = 0.5 #smaller number allows more collision
                    if dist < self.pt1.dist(self.pt2)*tolerance : 
                        self.isColliding = True

    def update(self) : 
        if self.isColliding : 
            self.del()
        elif self.time()==0 : 
            ICurve(self.pt1, self.pt2).clr(self.clr())
            
            if IRand.pct(95) : #branch 1
                dir = self.pt2.dif(self.pt1)
                angle = IRand.get(-PI/6, PI/6)
                dir.rot(angle)
                scale = IRand.get(0.90, 1.08)
                dir.mul(scale)
                MyAgent(self.pt2.cp(), dir).clr(self.clr())
            
            if IRand.pct(50) :  #branch 2
                dir = self.pt2.dif(self.pt1)
                angle = IRand.get(-PI/3*2, PI/3*2)
                dir.rot(angle)
                scale = IRand.get(0.90, 1.08)
                dir.mul(scale)
                MyAgent(self.pt2, dir).clr(self.clr())

The second code below puts polygon mesh boxes on the start points of the agents. Polygon mesh box is created by this method in IG class.
        IG.meshBox(corner, boxDir1, boxDir2, boxDir3)
The direction of the box is aligned to the direction of the agent and its size is randomized. The color is also randomized in gray scale.

add_library('igeo')

def setup() : 
    size(480, 360, IG.GL)
    IG.duration(80)
    MyAgent(IG.v(0, 0, 0), IG.v(1, 0, 0)).clr(0.5)


class MyAgent(IAgent) : 
    def __init__(self, pos, dir) : 
        self.pt1 = pos
        self.pt2 = pos.cp(dir)
        self.isColliding = False
    
    def interact(self, agents) : 
        for agent in agents : 
            if self.isColliding :
                return
            if isinstance(agent, MyAgent) : 
                if agent is not self : 
                    dist = agent.pt2.dist(self.pt2) #distance of end points
                    tolerance = 0.5 #smaller number allows more collision
                    if dist < self.pt1.dist(self.pt2)*tolerance : 
                        self.isColliding = True

    def update(self) : 
        if self.isColliding : 
            self.del()
        elif self.time()==0 : 
            boxScale = IRand.get(1.0, 2.5)
            boxDir1 = self.pt2.dif(self.pt1).mul(boxScale)
            boxDir2 = boxDir1.cp().rot(PI/2)
            boxDir3 = IG.v(0,0,boxDir1.len())
            corner = self.pt1.cp().sub(boxDir1.cp().div(2)).sub(boxDir2.cp().div(2))
            IG.meshBox(corner, boxDir1, boxDir2, boxDir3).clr(self.clr())
            
            if IRand.pct(95) : #branch 1
                dir = self.pt2.dif(self.pt1)
                angle = IRand.get(-PI/6, PI/6)
                dir.rot(angle)
                scale = IRand.get(0.90, 1.08)
                dir.mul(scale)
                gray = (self.red()+self.green()+self.blue())/3+ IRand.get(-0.05,0.05)
                MyAgent(self.pt2.cp(), dir).clr(gray)
            
            if IRand.pct(50) : #branch 2
                dir = self.pt2.dif(self.pt1)
                angle = IRand.get(-PI/3*2, PI/3*2)
                dir.rot(angle)
                scale = IRand.get(0.90, 1.08)
                dir.mul(scale)
                gray = (self.red()+self.green()+self.blue())/3 + IRand.get(-0.05,0.05)
                MyAgent(self.pt2, dir).clr(gray)

The next code manipulates a rotational axis of branching to create 3 dimensional branches. The agent class has a new data field axis as the rotational axis of branching rotation and also as one of edge directions of the box geometry.

add_library('igeo')

def setup() : 
    size(480, 360, IG.GL)
    IRand.init(2)
    IG.duration(40)
    MyAgent(IG.v(0, 0, 0), IG.v(1, 0, 0), IG.v(0, 0, 1)).clr(0.5)

class MyAgent(IAgent) : 
    def __init__(self, pos, dir, ax) : 
        self.pt1 = pos
        self.pt2 = pos.cp(dir)
        self.isColliding = False
        self.axis = ax
    
    def interact(self, agents) : 
        for agent in agents : 
            if self.isColliding :
                return
            if isinstance(agent, MyAgent) : 
                if agent is not self : 
                    dist = agent.pt2.dist(self.pt2) #distance of end points
                    tolerance = 0.5 #smaller number allows more collision
                    if dist < self.pt1.dist(self.pt2)*tolerance : 
                        self.isColliding = True
    
    def update(self) : 
        if self.isColliding : 
            self.del()
        elif self.time()==0 : 
            boxScale = IRand.get(1.0, 2.5)
            boxDir1 = self.pt2.dif(self.pt1).mul(boxScale)
            boxDir2 = boxDir1.cp().rot(self.axis, PI/2)
            boxDir3 = boxDir1.cross(boxDir2).len(boxDir1.len()) # perpendicular to boxDir1 and boxDir2
            corner = self.pt1.cp().sub(boxDir1.cp().div(2)).sub(boxDir2.cp().div(2))
            IG.meshBox(corner, boxDir1, boxDir2, boxDir3).clr(self.clr())
            self.axis = boxDir3
            
            if IRand.pct(95) : #branch 1
                dir = self.pt2.dif(self.pt1)
                angle = IRand.get(-PI/6, PI/6)
                dir.rot(self.axis, angle)
                scale = IRand.get(0.90, 1.08)
                dir.mul(scale)
                axis2 = self.axis.cp().rot(dir, IRand.get(-PI/6, PI/6))
                gray = (self.red()+self.green()+self.blue())/3+ IRand.get(-0.05,0.05)
                MyAgent(self.pt2.cp(), dir, axis2).clr(gray)
            
            if IRand.pct(50) : #branch 2
                dir = self.pt2.dif(self.pt1)
                angle = IRand.get(-PI/3*2, PI/3*2)
                dir.rot(self.axis, angle)
                scale = IRand.get(0.90, 1.08)
                dir.mul(scale)
                axis2 = self.axis.cp().rot(dir, IRand.get(-PI/6, PI/6))
                gray = (self.red()+self.green()+self.blue())/3 + IRand.get(-0.05,0.05)
                MyAgent(self.pt2, dir, axis2).clr(gray)

The code below takes a Rhino input file which contains one polygon mesh geometry and use this geometry as each agent's geometry instead of a box. The agent class has a new data field of mesh and the original template mesh geometry from the input file is passed to each agents. This template mesh geometry is copied, moved and rotated to the agent's orientation which is defined by pt1, pt2 and axis. The method to move and rotate is transform(xvec,yvec,zvec,translate).
        mesh.cp().transform(meshDir1, meshDir2, meshDir3, pt1);
The first three input arguments of IVec are new vectors to map the original X, Y and Z vectors onto. The fourth input argument of IVec is to move the whole geometry, just like adding a vector to the geometry.

The input file with a polygon mesh geometry used in the example is the following.

mesh1.3dm

add_library('igeo')

def setup() : 
    size(480, 360, IG.GL)
    IRand.init(2)
    IG.duration(40)
    IG.open("mesh1.3dm")
    mesh = IG.mesh(0)
    mesh.del() #hide the original
    MyAgent(mesh, IG.v(0, 0, 0), IG.v(1, 0, 0), IG.v(0, 0, 1)).clr(0.5)


class MyAgent(IAgent) : 
    def __init__(self, m, pos, dir, ax) : 
        self.pt1 = pos
        self.pt2 = pos.cp(dir)
        self.isColliding = False
        self.axis = ax
        self.mesh = m
    
    def interact(self, agents) : 
        for agent in agents : 
            if self.isColliding :
                return
            if isinstance(agent, MyAgent) : 
                if agent is not self : 
                    dist = agent.pt2.dist(self.pt2) #distance of end points
                    tolerance = 0.5 #smaller number allows more collision
                    if dist < self.pt1.dist(self.pt2)*tolerance : 
                        self.isColliding = True
    
    def update(self) : 
        if self.isColliding : 
            self.del()
        elif self.time()==0 : 
            meshScale = IRand.get(2.0, 5.0)
            meshDir1 = self.pt2.dif(self.pt1).mul(meshScale)
            meshDir2 = meshDir1.cp().rot(self.axis, PI/2)
            meshDir3 = meshDir1.cross(meshDir2).len(meshDir1.len()) # perpendicular to meshDir1 and meshDir2
            
            self.mesh.cp().transform(meshDir1, meshDir2, meshDir3, self.pt1).clr(self.clr())
            self.axis = meshDir3
            
            if IRand.pct(95) : #branch 1
                dir = self.pt2.dif(self.pt1)
                angle = IRand.get(-PI/3, PI/3)
                dir.rot(self.axis, angle)
                scale = IRand.get(0.90, 1.05)
                dir.mul(scale)
                axis2 = self.axis.cp().rot(dir, IRand.get(-PI/3, PI/3))
                gray = (self.red()+self.green()+self.blue())/3+ IRand.get(-0.05,0.05)
                MyAgent(self.mesh, self.pt2.cp(), dir, axis2).clr(gray)
            
            if IRand.pct(50) : #branch 2
                dir = self.pt2.dif(self.pt1)
                angle = IRand.get(-PI/3*2, PI/3*2)
                dir.rot(self.axis, angle)
                scale = IRand.get(0.90, 1.05)
                dir.mul(scale)
                axis2 = self.axis.cp().rot(dir, IRand.get(-PI/3, PI/3))
                gray = (self.red()+self.green()+self.blue())/3 + IRand.get(-0.05,0.05)
                MyAgent(self.mesh, self.pt2, dir, axis2).clr(gray)