The tutorial codes in this page extend the cell division and growth algorithm in the previous page by adding data structure of a triangular face which defines relationship of three cell agents in addition to a cell link which defines relationship of two, to explore surface formation through cell division and growth.

The first code in this page adds a new class CellFace which contains references to three cells and three links. Cell class also stores what CellFace instances contain the cell in the variable array field faces. CellLink also stores what CellFace instances contain the cell link in faces field.

The code below adds an algorithm to divide a face in Cell class's divide method. Conceptually the division of a face is done by dividing one of edge of the face triangle and insert a cell there and divide the original triangle into two with a new link in the middle as depicted in the diagram below. In the code, it's done by deleting the edge link to divide and the original face, and then add three new links and two new faces. When a cell contains only one link (only one line connecting two cells), it creates two links and one face with the new child cell to have a first face to be divided later.

add_library('igeo')

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    cell = Cell(IVec(0,0,0), 10)

class Cell(IParticle) : 
    growthDuration = 1400 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)
    
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(1.0,0,0)
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            child = self.createChild(IRand.dir()) # making a triangle loop
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        else : # strip state
            # divide one link
            dividingLink = self.links[IRand.getInt(0,len(self.links)-1)]
            if len(dividingLink.faces)==2 : 
                self.active = False
                return  # if link has two faces, skip division
            dir = dividingLink.oppositeDir(self) # vector to the other cell
            child = self.createChild(dir) # make a child toward the link direction
            c0 = dividingLink.oppositeCell(self) # opposite cell on the link
            
            f1 = dividingLink.faces[0] # existing face on the link
            c1 = f1.oppositeCell(dividingLink) # opposite cell from the link on the face
            # delete 1 link + 1 face, create 3 links + 2 faces
            CellLink(self, child)
            CellLink(c1,child)
            CellLink(c0,child)
            CellFace(self, c1, child)
            CellFace(child, c1, c0)
            dividingLink.delete()
            f1.delete()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
        self.line = ICurve(c1.pos(), c2.pos()).clr(1.0,0,0)
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.line.del() # delete line geometry
        self.del() # stop agent
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.cell1 = c1
        self.cell2 = c2
        self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
    
    def oppositeCell(self, l) : # find a cell opposite of a link edge in a triangle
        if self.cell1 is not l.cell1 and self.cell1 is not l.cell2 :
            return self.cell1
        if self.cell2 is not l.cell1 and self.cell2 is not l.cell2 :
            return self.cell2
        return self.cell3
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)

The code below has a sphere as the geometric representation of a cell instead of a circle. The code uses IG.meshSphere producing IMesh instance instead of using ISphere (NURBS surface sphere) to keep the code's memory usage lighter.

add_library('igeo')

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    cell = Cell(IVec(0,0,0), 10)
    cell.clr(0,0,1.0)

class Cell(IParticle) : 
    growthDuration = 1400 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.sphere = None
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)
    
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
        # update geometry
        if self.sphere is not None :
            self.sphere.del()
        self.sphere = IG.meshSphere(self.pos(), self.radius, 16).clr(self)
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            child = self.createChild(IRand.dir()) # making a triangle loop
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        else : # strip state
            # divide one link
            dividingLink = self.links[IRand.getInt(0,len(self.links)-1)]
            if len(dividingLink.faces)==2 : 
                self.active = False
                return  # if link has two faces, skip division
            dir = dividingLink.oppositeDir(self) # vector to the other cell
            child = self.createChild(dir) # make a child toward the link direction
            c0 = dividingLink.oppositeCell(self) # opposite cell on the link
            
            f1 = dividingLink.faces[0] # existing face on the link
            c1 = f1.oppositeCell(dividingLink) # opposite cell from the link on the face
            # delete 1 link + 1 face, create 3 links + 2 faces
            CellLink(self, child)
            CellLink(c1,child)
            CellLink(c0,child)
            CellFace(self, c1, child)
            CellFace(child, c1, c0)
            dividingLink.delete()
            f1.delete()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        child.hsb(self.hue()-0.05,1,0.8)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
        self.line = ICurve(c1.pos(), c2.pos()).clr(1.0,0,0)
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.line.del() # delete line geometry
        self.del() # stop agent
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.cell1 = c1
        self.cell2 = c2
        self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
    
    def oppositeCell(self, l) : # find a cell opposite of a link edge in a triangle
        if self.cell1 is not l.cell1 and self.cell1 is not l.cell2 :
            return self.cell1
        if self.cell2 is not l.cell1 and self.cell2 is not l.cell2 :
            return self.cell2
        return self.cell3
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)

The next code generates a polygon mesh out of cells, links and faces. A global variable mesh is defined and initialized in setup method, and then mesh faces are added at CellFace's constructor and also sometimes mesh faces are deleted at delete method. Cell class has a mesh vertex instance vertex for CellFace class to instantiate a mesh face (IFace) instance. Cell class and CellFace get nml method to calculate a normal vector of a face or of a vertex. To keep the normal direction consistent in adjacent faces, CellFace class adjust the order of cells using cellOrder method.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)

class Cell(IParticle) : 
    growthDuration = 1400 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)
    
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            child = self.createChild(IRand.dir()) # making a triangle loop
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        else : # strip state
            # divide one link
            dividingLink = self.links[IRand.getInt(0,len(self.links)-1)]
            if len(dividingLink.faces)==2 : 
                self.active = False
                return  # if link has two faces, skip division
            dir = dividingLink.oppositeDir(self) # vector to the other cell
            dir.projectToPlane(self.nml()) # division dir is projected on normal plane
            child = self.createChild(dir) # make a child toward the link direction
            c0 = dividingLink.oppositeCell(self) # opposite cell on the link
            
            f1 = dividingLink.faces[0] # existing face on the link
            c1 = f1.oppositeCell(dividingLink) # opposite cell from the link on the face
            # delete 1 link + 1 face, create 3 links + 2 faces
            CellLink(self, child)
            CellLink(c1,child)
            CellLink(c0,child)
            CellFace(self, c1, child)
            CellFace(child, c1, c0)
            dividingLink.delete()
            f1.delete()
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2
            self.cell2 = c1
            self.cell3 = c3
        else : 
            self.cell1 = c1
            self.cell2 = c2
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
    
    def oppositeCell(self, l) : # find a cell opposite of a link edge in a triangle
        if self.cell1 is not l.cell1 and self.cell1 is not l.cell2 :
            return self.cell1
        if self.cell2 is not l.cell1 and self.cell2 is not l.cell2 :
            return self.cell2
        return self.cell3
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

The previous face division algorithm always divides a peripheral edge which has only one adjacent face and never divides an interior edge which has two adjacent faces. The next code divides an interior edge in addition to peripheral edge. This is done in divide method and the use of the peripheral edge division or the interior edge division is determined by the number of faces the edge has (dividingLink.faces.size()).

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)

class Cell(IParticle) : 
    growthDuration = 1000 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)
    
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            child = self.createChild(IRand.dir()) # making a triangle loop
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        else : # strip state
            # divide one link
            dividingLink = self.links[IRand.getInt(0,len(self.links)-1)]
            
            dir = dividingLink.oppositeDir(self) # vector to the other cell
            dir.projectToPlane(self.nml()) # division dir is projected on normal plane
            child = self.createChild(dir) # make a child toward the link direction
            c0 = dividingLink.oppositeCell(self) # opposite cell on the link
            
            if len(dividingLink.faces)==1 : 
                f1 = dividingLink.faces[0] # existing face on the link
                c1 = f1.oppositeCell(dividingLink) # opposite cell from the link on the face
                # delete 1 link + 1 face, create 3 links + 2 faces
                CellLink(self, child)
                CellLink(c1,child)
                CellLink(c0,child)
                CellFace(self, c1, child)
                CellFace(child, c1, c0)
                dividingLink.delete()
                f1.delete()
            elif len(dividingLink.faces)==2 : 
                f1 = dividingLink.faces[0]
                f2 = dividingLink.faces[1]
                c1 = f1.oppositeCell(dividingLink)
                c2 = f2.oppositeCell(dividingLink)
                # delete 1 link + 1 face, create 4 links + 4 faces
                CellLink(self, child)
                CellLink(c1,child)
                CellLink(c2,child)
                CellLink(c0,child)
                CellFace(self, c1, child)
                CellFace(self, child, c2)
                CellFace(child, c1, c0)
                CellFace(child, c0, c2)
                dividingLink.delete()
                f1.delete()
                f2.delete()
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
    
    def oppositeCell(self, l) : # find a cell opposite of a link edge in a triangle
        if self.cell1 is not l.cell1 and self.cell1 is not l.cell2 :
            return self.cell1
        if self.cell2 is not l.cell1 and self.cell2 is not l.cell2 :
            return self.cell2
        return self.cell3
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

There is another way to divide a triangular face by insertion of one new node of a cell other than dividing one edge and dividing a triangle into two. If you insert a new cell node inside of a triangle and add three new edges, you can divide one face into three. The code below implements this face division algorithm.

The code generates closed polygon mesh by forming a tetrahedron after cells form a triangle. Once the network of cells, links and faces form a tetrahedron, the division algorithm never makes any holes and keep the whole mesh surface closed.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)

class Cell(IParticle) : 
    growthDuration = 1000 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)
    
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # divide one face
            dividingFace = self.faces[IRand.getInt(0,len(self.faces)-1)]
            c = dividingFace.oppositeCells(self)
            dir = dividingFace.center().dif(self.pos())
            dir.projectToPlane(self.nml()) 
            child = self.createChild(dir)
            # delete 1 face, create 3 links + 3 faces
            CellLink(self, child)
            CellLink(child, c[0])
            CellLink(child, c[1])
            CellFace(self, c[0], child)
            CellFace(self, child, c[1])
            CellFace(child, c[0], c[1])
            dividingFace.delete()
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCells(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeCell(self, l) : # find a cell opposite of a link edge in a triangle
        if self.cell1 is not l.cell1 and self.cell1 is not l.cell2 :
            return self.cell1
        if self.cell2 is not l.cell1 and self.cell2 is not l.cell2 :
            return self.cell2
        return self.cell3
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

The previous face division algorithm tends to generate sharp corners. The code below uses the interior edge division instead of the face division, which is the same algorithm with the one in the previous section.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)

class Cell(IParticle) : 
    growthDuration = 1000 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)
    
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # divide one link
            dividingLink = self.links[IRand.getInt(0,len(self.links)-1)]
            dir = dividingLink.oppositeDir(self)
            dir.projectToPlane(self.nml())
            child = self.createChild(dir)
            f1 = dividingLink.faces[0]
            f2 = dividingLink.faces[1]
            c0 = dividingLink.oppositeCell(self)
            c1 = f1.oppositeCell(dividingLink)
            c2 = f2.oppositeCell(dividingLink)
            # delete 1 link + 2 faces, create 4 links + 4 faces
            CellLink(self, child)
            CellLink(c0, child)
            CellLink(c1, child)
            CellLink(c2, child)
            CellFace(self, c1, child)
            CellFace(self, child, c2)
            CellFace(child, c1, c0)
            CellFace(child, c0, c2)
            dividingLink.delete()
            f1.delete()
            f2.delete()
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCells(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeCell(self, l) : # find a cell opposite of a link edge in a triangle
        if self.cell1 is not l.cell1 and self.cell1 is not l.cell2 :
            return self.cell1
        if self.cell2 is not l.cell1 and self.cell2 is not l.cell2 :
            return self.cell2
        return self.cell3
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

The face division and the edge division are not only way to divide a cell and recreate links and connection. To integrate different ways, the next codes implements an algorithm to select two edges on a dividing cell and insert a new child cell, three new edges and two new faces. The way of division is changed by the interval of two selected edges. If the interval is one, it's equivalent to the face division. If the interval is two, it's equivalent to the edge division.

This division algorithm is organized in several methods in addition to divide method. The two edges for division is selected at findDividingLinks method. To check the interval of two edges, the code needs to know which edge is next to each other and order of edges. For this purpose, sortLinks methods sorts the stored CellLink instances. insertChild method reconnect links and faces with a new child cell by deleting existing links and faces and creating new ones.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)

class Cell(IParticle) : 
    growthDuration = 1200 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)

        self.push(self.nml().len(20)) # pressure toward normal
  
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # pick two links and insert two faces between
            linkIdx = self.findDividingLinks()
            linksToChild = [] # links to reconnect to child
            for i in range(linkIdx[0], linkIdx[1]+1) : 
                linksToChild.append(self.links[i%len(self.links)])
            dir = IVec()
            for l in linksToChild :  # average of links1 dir
                dir.add(l.oppositeDir(self).unit())
            dir.projectToPlane(self.nml())
            child = self.createChild(dir)
            self.insertChild(child, linksToChild)
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child
    
    def sortLinks(self) : # sort links by following connected faces
        if len(self.links) <= 3 :
            return # no need to sort
        sorted = []
        currentLink = self.links[0]
        currentFace = currentLink.faces[0]
        for i in range(len(self.links)) : 
            sorted.append(currentLink)
            currentLink = currentFace.oppositeLink(currentLink.oppositeCell(self))
            currentFace = currentLink.oppositeFace(currentFace)
        self.links = sorted
    
    def findDividingLinks(self) :  # find splittable two links and return the index numbers
        linkInterval = len(self.links)//2 # interval between two links 
        
        self.sortLinks() # sort the order of links around cell
        idx = IRand.getInt(0,len(self.links)-1)
        return [ idx, idx+linkInterval ] # index can be larger than self.links.size()
    
    def insertChild(self, child, linksToChild) : # insert child cell and reconnect links
        removingFaces = []
        for f in self.faces : 
            for j in range(len(linksToChild)-1) : 
                if (f.contains(linksToChild[j]) and
                    f.contains(linksToChild[j+1])) : #a face between replacing links is removed
                    removingFaces.append(f)
                    break
        for i in range(1, len(linksToChild)-1) : 
            linksToChild[i].delete() # replacing links are once removed
            CellLink(linksToChild[i].oppositeCell(self), child) # then recreated
        cell1 = linksToChild[0].oppositeCell(self) # one of two dividing link cell
        cell2 = linksToChild[len(linksToChild)-1].oppositeCell(self) # another dividing link cell
        CellLink(child, cell1)
        CellLink(child, cell2)
        CellLink(self, child)
        for f in removingFaces : 
            f.delete() # replace face by deleting and recreating
            cells = f.oppositeCell(self)
            CellFace(child, cells[0], cells[1])
        CellFace(self, cell1, child)
        CellFace(self, child, cell2)

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeFace(self, f) : # find other face on the link
        if len(self.faces)!=2 :
            return None
        if self.faces[0] is f :
            return self.faces[1]
        return self.faces[0]
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCell(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeLink(self, c) : # find opposite link of a cell in a triangle
        if not self.link1.contains(c) :
            return self.link1
        if not self.link2.contains(c) :
            return self.link2
        if not self.link3.contains(c) :
            return self.link3
        return None
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

The next code brings a rule to limit number of links in a cell. The number of links is checked in findDividingLinks method not only with the dividing cell but also with adjacent cells which would increase links if new links and faces are inserted. findDividingLinks also searches a pair of edges which satisfy the maximum link number condition by checking all possible pairs with the specified link interval until it finds one.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)

class Cell(IParticle) : 
    growthDuration = 20000 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    maxLink = 6 #limit on number of links per cell
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)

        self.push(self.nml().len(20)) # pressure toward normal
  
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # pick two links and insert two faces between
            linkIdx = self.findDividingLinks()
            if linkIdx is None : # no valid edge to split. skip division.
                self.active = False
                return 
            linksToChild = [] # links to reconnect to child
            for i in range(linkIdx[0], linkIdx[1]+1) : 
                linksToChild.append(self.links[i%len(self.links)])
            dir = IVec()
            for l in linksToChild :  # average of links1 dir
                dir.add(l.oppositeDir(self).unit())
            dir.projectToPlane(self.nml())
            child = self.createChild(dir)
            self.insertChild(child, linksToChild)
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child
    
    def sortLinks(self) : # sort links by following connected faces
        if len(self.links) <= 3 :
            return # no need to sort
        sorted = []
        currentLink = self.links[0]
        currentFace = currentLink.faces[0]
        for i in range(len(self.links)) : 
            sorted.append(currentLink)
            currentLink = currentFace.oppositeLink(currentLink.oppositeCell(self))
            currentFace = currentLink.oppositeFace(currentFace)
        self.links = sorted
    
    def findDividingLinks(self) :  # find splittable two links and return the index numbers
        linkInterval = len(self.links)//2 # interval between two links 
        
        self.sortLinks() # sort the order of links around cell
        if linkInterval==1 and len(self.links) >= Cell.maxLink  : # dividing links next to each other adds one more link

            return None
        idx = IRand.getInt(0,len(self.links)-1)
        for i in range(len(self.links)) : # check all pairs
            c1 = self.links[(i+idx)%len(self.links)].oppositeCell(self)
            c2 = self.links[(i+idx+linkInterval)%len(self.links)].oppositeCell(self) 
            if len(c1.links) < Cell.maxLink and len(c2.links) < Cell.maxLink : # division adds one link on the end of each dividing link
                return [ i+idx, i+idx+linkInterval ] # index can be larger than self.links.size()
        return None
    
    def insertChild(self, child, linksToChild) : # insert child cell and reconnect links
        removingFaces = []
        for f in self.faces : 
            for j in range(len(linksToChild)-1) : 
                if (f.contains(linksToChild[j]) and
                    f.contains(linksToChild[j+1])) : #a face between replacing links is removed
                    removingFaces.append(f)
                    break
        for i in range(1, len(linksToChild)-1) : 
            linksToChild[i].delete() # replacing links are once removed
            CellLink(linksToChild[i].oppositeCell(self), child) # then recreated
        cell1 = linksToChild[0].oppositeCell(self) # one of two dividing link cell
        cell2 = linksToChild[len(linksToChild)-1].oppositeCell(self) # another dividing link cell
        CellLink(child, cell1)
        CellLink(child, cell2)
        CellLink(self, child)
        for f in removingFaces : 
            f.delete() # replace face by deleting and recreating
            cells = f.oppositeCell(self)
            CellFace(child, cells[0], cells[1])
        CellFace(self, cell1, child)
        CellFace(self, child, cell2)

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeFace(self, f) : # find other face on the link
        if len(self.faces)!=2 :
            return None
        if self.faces[0] is f :
            return self.faces[1]
        return self.faces[0]
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCell(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeLink(self, c) : # find opposite link of a cell in a triangle
        if not self.link1.contains(c) :
            return self.link1
        if not self.link2.contains(c) :
            return self.link2
        if not self.link3.contains(c) :
            return self.link3
        return None
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

Although the limitation of 6 maximum links tends to bring a certain regularity in the surface formation, it also tends to limit the growth because number of links can be easily maxed out and a dividable pair of edges cannot be found. To ease the condition, the next code accepts one more link randomly with a certain probability (extraLinkAllowance field). This change happens in findDividingLinks method.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)
    cell.active = True

class Cell(IParticle) : 
    growthDuration = 2000 #duration of growth and division
    growthInterval = 10
    divisionInterval = 50
    maxRadius = 100
    growthSpeed=0.1
    maxLink = 6 #limit on number of links per cell
    extraLinkAllowance = 20 #percentage to accept one more link than maxLink
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)

        self.push(self.nml().len(20)) # pressure toward normal
  
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
                if IRand.pct(50) : # random division
                    self.active = True
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # pick two links and insert two faces between
            linkIdx = self.findDividingLinks()
            if linkIdx is None : # no valid edge to split. skip division.
                self.active = False
                return 
            linksToChild = [] # links to reconnect to child
            for i in range(linkIdx[0], linkIdx[1]+1) : 
                linksToChild.append(self.links[i%len(self.links)])
            dir = IVec()
            for l in linksToChild :  # average of links1 dir
                dir.add(l.oppositeDir(self).unit())
            dir.projectToPlane(self.nml())
            child = self.createChild(dir)
            self.insertChild(child, linksToChild)
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        return child
    
    def sortLinks(self) : # sort links by following connected faces
        if len(self.links) <= 3 :
            return # no need to sort
        sorted = []
        currentLink = self.links[0]
        currentFace = currentLink.faces[0]
        for i in range(len(self.links)) : 
            sorted.append(currentLink)
            currentLink = currentFace.oppositeLink(currentLink.oppositeCell(self))
            currentFace = currentLink.oppositeFace(currentFace)
        self.links = sorted
    
    def findDividingLinks(self) :  # find splittable two links and return the index numbers
        linkInterval = len(self.links)//2 # interval between two links 
        
        self.sortLinks() # sort the order of links around cell
        allowExtra = IRand.pct(Cell.extraLinkAllowance)
        if (linkInterval==1 and # dividing links next to each other adds one more link
            not (len(self.links) < Cell.maxLink or allowExtra and len(self.links) < Cell.maxLink+1) ) : 
            return None
        idx = IRand.getInt(0,len(self.links)-1)
        for i in range(len(self.links)) : # check all pairs
            c1 = self.links[(i+idx)%len(self.links)].oppositeCell(self)
            c2 = self.links[(i+idx+linkInterval)%len(self.links)].oppositeCell(self) 
            if (len(c1.links) < Cell.maxLink and len(c2.links) < Cell.maxLink or # division adds one link on the end of each dividing link
                allowExtra and len(c1.links) < Cell.maxLink+1 and len(c2.links) < Cell.maxLink+1 ) : 
                return [ i+idx, i+idx+linkInterval ] # index can be larger than self.links.size()
        return None
    
    def insertChild(self, child, linksToChild) : # insert child cell and reconnect links
        removingFaces = []
        for f in self.faces : 
            for j in range(len(linksToChild)-1) : 
                if (f.contains(linksToChild[j]) and
                    f.contains(linksToChild[j+1])) : #a face between replacing links is removed
                    removingFaces.append(f)
                    break
        for i in range(1, len(linksToChild)-1) : 
            linksToChild[i].delete() # replacing links are once removed
            CellLink(linksToChild[i].oppositeCell(self), child) # then recreated
        cell1 = linksToChild[0].oppositeCell(self) # one of two dividing link cell
        cell2 = linksToChild[len(linksToChild)-1].oppositeCell(self) # another dividing link cell
        CellLink(child, cell1)
        CellLink(child, cell2)
        CellLink(self, child)
        for f in removingFaces : 
            f.delete() # replace face by deleting and recreating
            cells = f.oppositeCell(self)
            CellFace(child, cells[0], cells[1])
        CellFace(self, cell1, child)
        CellFace(self, child, cell2)

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeFace(self, f) : # find other face on the link
        if len(self.faces)!=2 :
            return None
        if self.faces[0] is f :
            return self.faces[1]
        return self.faces[0]
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCell(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeLink(self, c) : # find opposite link of a cell in a triangle
        if not self.link1.contains(c) :
            return self.link1
        if not self.link2.contains(c) :
            return self.link2
        if not self.link3.contains(c) :
            return self.link3
        return None
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

The next three codes below bring a type of division control by activation parameter (active field in the Cell class) and controlling it by a parent cell when it creates a new child cell in division method. The next code use the cascading logic to activate a child cell and deactivate the parent cell after it finishes division, in the same way with the example in the previous page.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)
    cell.active = True

class Cell(IParticle) : 
    growthDuration = 10000 #duration of growth and division
    growthInterval = 10
    divisionInterval = 20
    maxRadius = 100
    growthSpeed=0.1
    maxLink = 6 #limit on number of links per cell
    extraLinkAllowance = 0 #percentage to accept one more link than maxLink
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(100) # constant force
            self.push(dif)

        self.push(self.nml().len(20)) # pressure toward normal
  
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # pick two links and insert two faces between
            linkIdx = self.findDividingLinks()
            if linkIdx is None : # no valid edge to split. skip division.
                self.active = False
                return 
            linksToChild = [] # links to reconnect to child
            for i in range(linkIdx[0], linkIdx[1]+1) : 
                linksToChild.append(self.links[i%len(self.links)])
            dir = IVec()
            for l in linksToChild :  # average of links1 dir
                dir.add(l.oppositeDir(self).unit())
            dir.projectToPlane(self.nml())
            child = self.createChild(dir)
            self.insertChild(child, linksToChild)
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        self.active = False #reset activation
        child.active = True #activate child always
        return child
    
    def sortLinks(self) : # sort links by following connected faces
        if len(self.links) <= 3 :
            return # no need to sort
        sorted = []
        currentLink = self.links[0]
        currentFace = currentLink.faces[0]
        for i in range(len(self.links)) : 
            sorted.append(currentLink)
            currentLink = currentFace.oppositeLink(currentLink.oppositeCell(self))
            currentFace = currentLink.oppositeFace(currentFace)
        self.links = sorted
    
    def findDividingLinks(self) :  # find splittable two links and return the index numbers
        linkInterval = 1 # face division only 
        
        self.sortLinks() # sort the order of links around cell
        allowExtra = IRand.pct(Cell.extraLinkAllowance)
        if (linkInterval==1 and # dividing links next to each other adds one more link
            not (len(self.links) < Cell.maxLink or allowExtra and len(self.links) < Cell.maxLink+1) ) : 
            return None
        idx = IRand.getInt(0,len(self.links)-1)
        for i in range(len(self.links)) : # check all pairs
            c1 = self.links[(i+idx)%len(self.links)].oppositeCell(self)
            c2 = self.links[(i+idx+linkInterval)%len(self.links)].oppositeCell(self) 
            if (len(c1.links) < Cell.maxLink and len(c2.links) < Cell.maxLink or # division adds one link on the end of each dividing link
                allowExtra and len(c1.links) < Cell.maxLink+1 and len(c2.links) < Cell.maxLink+1 ) : 
                return [ i+idx, i+idx+linkInterval ] # index can be larger than self.links.size()
        return None
    
    def insertChild(self, child, linksToChild) : # insert child cell and reconnect links
        removingFaces = []
        for f in self.faces : 
            for j in range(len(linksToChild)-1) : 
                if (f.contains(linksToChild[j]) and
                    f.contains(linksToChild[j+1])) : #a face between replacing links is removed
                    removingFaces.append(f)
                    break
        for i in range(1, len(linksToChild)-1) : 
            linksToChild[i].delete() # replacing links are once removed
            CellLink(linksToChild[i].oppositeCell(self), child) # then recreated
        cell1 = linksToChild[0].oppositeCell(self) # one of two dividing link cell
        cell2 = linksToChild[len(linksToChild)-1].oppositeCell(self) # another dividing link cell
        CellLink(child, cell1)
        CellLink(child, cell2)
        CellLink(self, child)
        for f in removingFaces : 
            f.delete() # replace face by deleting and recreating
            cells = f.oppositeCell(self)
            CellFace(child, cells[0], cells[1])
        CellFace(self, cell1, child)
        CellFace(self, child, cell2)

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeFace(self, f) : # find other face on the link
        if len(self.faces)!=2 :
            return None
        if self.faces[0] is f :
            return self.faces[1]
        return self.faces[0]
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCell(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeLink(self, c) : # find opposite link of a cell in a triangle
        if not self.link1.contains(c) :
            return self.link1
        if not self.link2.contains(c) :
            return self.link2
        if not self.link3.contains(c) :
            return self.link3
        return None
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

The division algorithm in the next code activates a child cell and also keeps the parent cell activated to cause branching, like the example in the previous page. Parent cells are maintained to be active only in 10% probability (90% deactivated) in the code not to generate exponentially many branches.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)
    cell.active = True

class Cell(IParticle) : 
    growthDuration = 2000 #duration of growth and division
    growthInterval = 10
    divisionInterval = 20
    maxRadius = 100
    growthSpeed=0.1
    maxLink = 8 #limit on number of links per cell
    extraLinkAllowance = 0 #percentage to accept one more link than maxLink
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(100) # constant force
            self.push(dif)

        self.push(self.nml().len(20)) # pressure toward normal
  
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # pick two links and insert two faces between
            linkIdx = self.findDividingLinks()
            if linkIdx is None : # no valid edge to split. skip division.
                self.active = False
                return 
            linksToChild = [] # links to reconnect to child
            for i in range(linkIdx[0], linkIdx[1]+1) : 
                linksToChild.append(self.links[i%len(self.links)])
            dir = IVec()
            for l in linksToChild :  # average of links1 dir
                dir.add(l.oppositeDir(self).unit())
            dir.projectToPlane(self.nml())
            child = self.createChild(dir)
            self.insertChild(child, linksToChild)
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        if IRand.pct(90) :
            self.active = False #10% stay active
        child.active = True #activate child always
        return child
    
    def sortLinks(self) : # sort links by following connected faces
        if len(self.links) <= 3 :
            return # no need to sort
        sorted = []
        currentLink = self.links[0]
        currentFace = currentLink.faces[0]
        for i in range(len(self.links)) : 
            sorted.append(currentLink)
            currentLink = currentFace.oppositeLink(currentLink.oppositeCell(self))
            currentFace = currentLink.oppositeFace(currentFace)
        self.links = sorted
    
    def findDividingLinks(self) :  # find splittable two links and return the index numbers
        linkInterval = 1 # interaval between two links 
        
        self.sortLinks() # sort the order of links around cell
        allowExtra = IRand.pct(Cell.extraLinkAllowance)
        if (linkInterval==1 and # dividing links next to each other adds one more link
            not (len(self.links) < Cell.maxLink or allowExtra and len(self.links) < Cell.maxLink+1) ) : 
            return None
        idx = IRand.getInt(0,len(self.links)-1)
        for i in range(len(self.links)) : # check all pairs
            c1 = self.links[(i+idx)%len(self.links)].oppositeCell(self)
            c2 = self.links[(i+idx+linkInterval)%len(self.links)].oppositeCell(self) 
            if (len(c1.links) < Cell.maxLink and len(c2.links) < Cell.maxLink or # division adds one link on the end of each dividing link
                allowExtra and len(c1.links) < Cell.maxLink+1 and len(c2.links) < Cell.maxLink+1 ) : 
                return [ i+idx, i+idx+linkInterval ] # index can be larger than self.links.size()
        return None
    
    def insertChild(self, child, linksToChild) : # insert child cell and reconnect links
        removingFaces = []
        for f in self.faces : 
            for j in range(len(linksToChild)-1) : 
                if (f.contains(linksToChild[j]) and
                    f.contains(linksToChild[j+1])) : #a face between replacing links is removed
                    removingFaces.append(f)
                    break
        for i in range(1, len(linksToChild)-1) : 
            linksToChild[i].delete() # replacing links are once removed
            CellLink(linksToChild[i].oppositeCell(self), child) # then recreated
        cell1 = linksToChild[0].oppositeCell(self) # one of two dividing link cell
        cell2 = linksToChild[len(linksToChild)-1].oppositeCell(self) # another dividing link cell
        CellLink(child, cell1)
        CellLink(child, cell2)
        CellLink(self, child)
        for f in removingFaces : 
            f.delete() # replace face by deleting and recreating
            cells = f.oppositeCell(self)
            CellFace(child, cells[0], cells[1])
        CellFace(self, cell1, child)
        CellFace(self, child, cell2)

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeFace(self, f) : # find other face on the link
        if len(self.faces)!=2 :
            return None
        if self.faces[0] is f :
            return self.faces[1]
        return self.faces[0]
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCell(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeLink(self, c) : # find opposite link of a cell in a triangle
        if not self.link1.contains(c) :
            return self.link1
        if not self.link2.contains(c) :
            return self.link2
        if not self.link3.contains(c) :
            return self.link3
        return None
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

The code below introduces two parameters in Cell class; stayActive and activateChild. They control whether a parent keeps activated and whether a child cell is activated when it's created in cell division. Now it's possible to control the activation control type to be cascading, branching or others outside of createChild method by setting the parameters in a cell.

add_library('igeo')

mesh = None # global variable of mesh geometry

def setup() : 
    size(480,360,IG.GL)
    IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
    mesh = IMesh().clr(0.7,0,0)
    cell = Cell(IVec(0,0,0), 10)
    cell.stayActive = True
    cell.activateChild = True
    cell.active = True

class Cell(IParticle) : 
    growthDuration = 1800 #duration of growth and division
    growthInterval = 10
    divisionInterval = 20
    maxRadius = 100
    growthSpeed=0.1
    maxLink = 6 #limit on number of links per cell
    extraLinkAllowance = 10 #percentage to accept one more link than maxLink
    
    def __init__(self, pos, rad) :
        IParticle.__init__(self, pos, IVec(0,0,0))
        self.radius = rad
        self.links = [] #store links
        self.faces = [] #store faces
        self.active = False
        self.stayActive = False #activation control variables
        self.activateChild = False
        self.vertex = IVertex(self.pos())
        self.fric(0.2)
  
    def interact(self, agents) : 
        neighborCenter = IVec(0,0,0)
        neighborCount=0
        neighborDist = self.radius*4 
        for a in agents : 
            if isinstance(a, Cell) : 
                if a is not self : 
                    # push if closer than two radii
                    if a.pos().dist(self.pos()) < a.radius+self.radius : 
                        dif = a.pos().dif(self.pos())
                        dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
                        a.push(dif)
                    # count neighbors and calculate their center
                    if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
                        neighborCenter.add(a.pos())
                        neighborCount+=1
        if neighborCount >= 1 : # push from center of neighbors
            neighborCenter.div(neighborCount)
            dif = self.pos().dif(neighborCenter).len(50) # constant force
            self.push(dif)

        self.push(self.nml().len(100)) # pressure toward normal
  
    def update(self) : 
        if IG.time() < Cell.growthDuration : 
            if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
                if self.active : # divide when active flag is on
                    self.divide()
            if self.time()%Cell.growthInterval==0 : 
                self.grow()
            if self.active : # update color by being active
                self.clr(1.0,0.5,0)
            else :
                self.clr(mesh.clr())
        self.vertex.pos().set(self.pos()) # update mesh vertex position
        self.vertex.nml(self.nml()) # update mesh vertex normal
    
    def grow(self) : # growing cell size
        if self.radius < Cell.maxRadius : 
            self.radius += Cell.growthSpeed
  
    def divide(self) : # cell division
        if len(self.links)==0 : # dot state
            child = self.createChild(IRand.dir())
            CellLink(self, child) # add one link
        elif len(self.links)==1 : # line state 
            # making a triangle loop
            child = self.createChild(IRand.dir())
            CellLink(child, self.links[0].cell1)
            CellLink(child, self.links[0].cell2)
            CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
        elif len(self.links)==2 : # triangle
            # making a tetrahedron enclosure
            child = self.createChild(IRand.dir())
            f = self.faces[0]
            center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
            if f.center().dif(center).dot(f.nml())<0 :
                f.flipNml() # adjust normal to be outward
            CellLink(f.cell1, child)
            CellLink(f.cell2, child)
            CellLink(f.cell3, child)
            f1 = CellFace(f.cell1, f.cell2, child)
            f2 = CellFace(f.cell2, f.cell3, child)
            f3 = CellFace(f.cell3, f.cell1, child)
            if f1.center().dif(center).dot(f1.nml())<0 : 
                f1.flipNml()
            if f2.center().dif(center).dot(f2.nml())<0 : 
                f2.flipNml()
            if f3.center().dif(center).dot(f3.nml())<0 : 
                f3.flipNml()
        else : # link num > 3
            # pick two links and insert two faces between
            linkIdx = self.findDividingLinks()
            if linkIdx is None : # no valid edge to split. skip division.
                self.active = False
                return 
            linksToChild = [] # links to reconnect to child
            for i in range(linkIdx[0], linkIdx[1]+1) : 
                linksToChild.append(self.links[i%len(self.links)])
            dir = IVec()
            for l in linksToChild :  # average of links1 dir
                dir.add(l.oppositeDir(self).unit())
            dir.projectToPlane(self.nml())
            child = self.createChild(dir)
            self.insertChild(child, linksToChild)
    
    def nml(self) : # calc vertex normal from face normal
        if len(self.faces)==0 :
            return IVec(0,0,1)
        n = IVec(0,0,0)
        for f in self.faces : 
            n.add(f.nml())
        return n.unit()
    
    def createChild(self, dir) : 
        self.radius *= 0.5 #make both cell size half
        dir.len(self.radius)
        child = Cell(self.pos().cp(dir), self.radius)
        self.pos().sub(dir) # move to the other way
        if self.stayActive : #activation control
            self.active = True
        else :
            self.active = False
        if self.activateChild :
            child.active = True
        else :
            child.active = False
        child.stayActive = self.stayActive #activation behavior parameters are passed to child
        child.activateChild = self.activateChild
        return child
  
    def sortLinks(self) : # sort links by following connected faces
        if len(self.links) <= 3 :
            return # no need to sort
        sorted = []
        currentLink = self.links[0]
        currentFace = currentLink.faces[0]
        for i in range(len(self.links)) : 
            sorted.append(currentLink)
            currentLink = currentFace.oppositeLink(currentLink.oppositeCell(self))
            currentFace = currentLink.oppositeFace(currentFace)
        self.links = sorted
    
    def findDividingLinks(self) :  # find splittable two links and return the index numbers
        linkInterval = len(self.links)//2 # interaval between two links 
        
        self.sortLinks() # sort the order of links around cell
        allowExtra = IRand.pct(Cell.extraLinkAllowance)
        if (linkInterval==1 and # dividing links next to each other adds one more link
            not (len(self.links) < Cell.maxLink or allowExtra and len(self.links) < Cell.maxLink+1) ) : 
            return None
        idx = IRand.getInt(0,len(self.links)-1)
        for i in range(len(self.links)) : # check all pairs
            c1 = self.links[(i+idx)%len(self.links)].oppositeCell(self)
            c2 = self.links[(i+idx+linkInterval)%len(self.links)].oppositeCell(self) 
            if (len(c1.links) < Cell.maxLink and len(c2.links) < Cell.maxLink or # division adds one link on the end of each dividing link
                allowExtra and len(c1.links) < Cell.maxLink+1 and len(c2.links) < Cell.maxLink+1 ) : 
                return [ i+idx, i+idx+linkInterval ] # index can be larger than self.links.size()
        return None
    
    def insertChild(self, child, linksToChild) : # insert child cell and reconnect links
        removingFaces = []
        for f in self.faces : 
            for j in range(len(linksToChild)-1) : 
                if (f.contains(linksToChild[j]) and
                    f.contains(linksToChild[j+1])) : #a face between replacing links is removed
                    removingFaces.append(f)
                    break
        for i in range(1, len(linksToChild)-1) : 
            linksToChild[i].delete() # replacing links are once removed
            CellLink(linksToChild[i].oppositeCell(self), child) # then recreated
        cell1 = linksToChild[0].oppositeCell(self) # one of two dividing link cell
        cell2 = linksToChild[len(linksToChild)-1].oppositeCell(self) # another dividing link cell
        CellLink(child, cell1)
        CellLink(child, cell2)
        CellLink(self, child)
        for f in removingFaces : 
            f.delete() # replace face by deleting and recreating
            cells = f.oppositeCell(self)
            CellFace(child, cells[0], cells[1])
        CellFace(self, cell1, child)
        CellFace(self, child, cell2)

class CellLink(IAgent) : # a link to connect 2 cells with spring force
    maxForce = 100
    
    def __init__(self, c1, c2) : 
        IAgent.__init__(self)
        self.cell1 = c1 
        self.cell2 = c2
        self.cell1.links.append(self) # register this link to cells{
        self.cell2.links.append(self)
        self.faces = []
    
    def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

    def contains(self, c) : #check if link contains the cell
        if c is self.cell1 or c is self.cell2 : 
            return True
        return False
    
    def delete(self) : 
        self.cell1.links.remove(self) # unregister from cells
        self.cell2.links.remove(self)
        self.del() # stop agent
  
    def oppositeFace(self, f) : # find other face on the link
        if len(self.faces)!=2 :
            return None
        if self.faces[0] is f :
            return self.faces[1]
        return self.faces[0]
  
    def oppositeCell(self, c) : # find other cell on the link
        if self.cell1 is c : 
            return self.cell2
        if self.cell2 is c :
            return self.cell1
        print("Link does not contain the input cell")
        return None 
  
    def oppositeDir(self, c) : #calculate a vector to the other cell
        return self.oppositeCell(c).pos().dif(c.pos())

class CellFace : # a triangle grouping 3 cells and 3 links
    
    def __init__(self, c1, c2, c3) : 
        self.link1 = self.findLink(c1, c2)
        # keep order of cells in right hand screw for consistent normal
        if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
            self.cell1 = c2 
            self.cell2 = c1 
            self.cell3 = c3
        else : 
            self.cell1 = c1 
            self.cell2 = c2 
            self.cell3 = c3
        self.link1 = self.findLink(self.cell1, self.cell2)
        self.link2 = self.findLink(self.cell2, self.cell3)
        self.link3 = self.findLink(self.cell3, self.cell1)
        self.cell1.faces.append(self) # register this face to cells and links
        self.cell2.faces.append(self)
        self.cell3.faces.append(self)
        self.link1.faces.append(self)
        self.link2.faces.append(self)
        self.link3.faces.append(self)
        self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
        mesh.addFace(self.face) # add triangular face to the mesh
    
    def center(self) : # calc center of triangle
        return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
  
    def nml(self) : # normal vector
        return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
    
    def flipNml(self) : # flip normal
        tmp = self.cell2
        tmpLink = self.link2
        self.cell2 = self.cell3
        self.link2 = self.link3
        self.cell3 = tmp
        self.link3 = tmpLink

    # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
    def cellOrder(self, c1, c2) : 
        if ( self.cell1 is c1 and self.cell2 is c2 or 
             self.cell2 is c1 and self.cell3 is c2 or 
             self.cell3 is c1 and self.cell1 is c2 ) : 
                return True
        return False
    
    def contains(self, link) : # check if the link is contained
        return self.link1 is link or self.link2 is link or self.link3 is link
    
    def findLink(self, c1, c2) : # find a link between 2 cells
        for l in c1.links : 
            if l.contains(c2) : 
                return l
        print("link not found")
        return None
  
    def oppositeCell(self, c) : # find 2 cells opposite of a cell in a triangle
        if self.cell1 is c :
            return  [ self.cell2, self.cell3 ]
        if self.cell2 is c :
            return  [ self.cell3, self.cell1 ]
        if self.cell3 is c :
            return  [ self.cell1, self.cell2 ]
        return None
  
    def oppositeLink(self, c) : # find opposite link of a cell in a triangle
        if not self.link1.contains(c) :
            return self.link1
        if not self.link2.contains(c) :
            return self.link2
        if not self.link3.contains(c) :
            return self.link3
        return None
    
    def delete(self) : 
        self.cell1.faces.remove(self) # unregister self from cells and links
        self.cell2.faces.remove(self)
        self.cell3.faces.remove(self)
        self.link1.faces.remove(self)
        self.link2.faces.remove(self)
        self.link3.faces.remove(self)
        mesh.deleteFace(self.face)

The code below import a file which contains several meshes and build a network of cells, cell links and cell faces for each mesh. The method cellsFromMesh generates cells. You can specify which cells to be fixed adn which cells to be activated by locating points in spefici layers. The codes of Cell, CellLink, CellFace class add some fields and methods to deal with exsisting mesh and both open and closed mesh.

The sample mesh file used in the example is below.

mesh2.3dm

Warning:  include(codepy/igeo_tutorial56_13/igeo_tutorial56_13.html): failed to open stream: No such file or directory in /home/mj7789dybiu5/public_html/igeo/tutorial/tutorial.php on line 66



Warning:  include(): Failed opening 'codepy/igeo_tutorial56_13/igeo_tutorial56_13.html' for inclusion (include_path='.:/opt/alt/php56/usr/share/pear:/opt/alt/php56/usr/share/php') in /home/mj7789dybiu5/public_html/igeo/tutorial/tutorial.php on line 66