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.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  Cell cell = new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 1400; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(1.0,0,0); }
    }
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1); // making a triangle loop
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else{ // strip state
       // divide one link
      CellLink dividingLink = links.get(IRand.getInt(0,links.size()-1));
      if(dividingLink.faces.size()==2){
        active = false;
        return; // if link has two faces, skip division
      }
      IVec dir = dividingLink.oppositeDir(this); // vector to the other cell
      Cell child = createChild(dir); // make a child toward the link direction
      Cell c0 = dividingLink.oppositeCell(this); // opposite cell on the link
      
      CellFace f1 = dividingLink.faces.get(0); // existing face on the link
      Cell c1 = f1.oppositeCell(dividingLink); // opposite cell from the link on the face
      // delete 1 link + 1 face, create 3 links + 2 faces
      new CellLink(this, child);
      new CellLink(c0,child);
      new CellLink(c1,child);
      new CellFace(this, c1, child);
      new CellFace(child, c1, c0);
      dividingLink.del();
      f1.del();
    }
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  ICurve line;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
    line = new ICurve(c1.pos(), c2.pos()).clr(1.,0,0);
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    line.del(); // delete line geometry
    super.del(); // stop agent
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    cell1 = c1;
    cell2 = c2;
    cell3 = c3;
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell oppositeCell(CellLink l){ // find a cell opposite of a link edge in a triangle
    if(cell1!=l.cell1 && cell1!=l.cell2){ return cell1; }
    if(cell2!=l.cell1 && cell2!=l.cell2){ return cell2; }
    return cell3;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
  }
}

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.

import igeo.*;
import processing.opengl.*;

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  Cell cell = new Cell(new IVec(0,0,0), 10);
  cell.clr(0,0,1.0);
}

class Cell extends IParticle{
  int growthDuration = 1400; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IMesh sphere;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
    }
    // update geometry
    if(sphere!=null) sphere.del();
    sphere = IG.meshSphere(pos(), radius, 16).clr(clr());
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1); // making a triangle loop
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else{ // strip state
       // divide one link
      CellLink dividingLink = links.get(IRand.getInt(0,links.size()-1));
      if(dividingLink.faces.size()==2){
        active = false;
        return;  // if link has two faces, skip division
      }
      IVec dir = dividingLink.oppositeDir(this); // vector to the other cell
      Cell child = createChild(dir); // make a child toward the link direction
      Cell c0 = dividingLink.oppositeCell(this); // opposite cell on the link
      
      CellFace f1 = dividingLink.faces.get(0); // existing face on the link
      Cell c1 = f1.oppositeCell(dividingLink); // opposite cell from the link on the face
      // delete 1 link + 1 face, create 3 links + 2 faces
      new CellLink(this, child);
      new CellLink(c1,child);
      new CellLink(c0,child);
      new CellFace(this, c1, child);
      new CellFace(child, c1, c0);
      dividingLink.del();
      f1.del();
    }
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    child.hsb(hue()-0.05,1,0.8);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  ICurve line;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
    line = new ICurve(c1.pos(), c2.pos()).clr(1.,0,0);
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    line.del(); // delete line geometry
    super.del(); // stop agent
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    cell1 = c1;
    cell2 = c2;
    cell3 = c3;
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell oppositeCell(CellLink l){ // find a cell opposite of a link edge in a triangle
    if(cell1!=l.cell1 && cell1!=l.cell2){ return cell1; }
    if(cell2!=l.cell1 && cell2!=l.cell2){ return cell2; }
    return cell3;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
  }
}

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

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 1400; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1); // making a triangle loop
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else{ // strip state
       // divide one link
      CellLink dividingLink = links.get(IRand.getInt(0,links.size()-1));
      if(dividingLink.faces.size()==2){
        active = false;
        return;  // if link has two faces, skip division
      }
      IVec dir = dividingLink.oppositeDir(this); // vector to the other cell
      dir.projectToPlane(nml()); // division dir is projected on normal plane
      Cell child = createChild(dir); // make a child toward the link direction
      Cell c0 = dividingLink.oppositeCell(this); // opposite cell on the link
      
      CellFace f1 = dividingLink.faces.get(0); // existing face on the link
      Cell c1 = f1.oppositeCell(dividingLink); // opposite cell from the link on the face
      // delete 1 link + 1 face, create 3 links + 2 faces
      new CellLink(this, child);
      new CellLink(c1,child);
      new CellLink(c0,child);
      new CellFace(this, c1, child);
      new CellFace(child, c1, c0);
      dividingLink.del();
      f1.del();
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell oppositeCell(CellLink l){ // find a cell opposite of a link edge in a triangle
    if(cell1!=l.cell1 && cell1!=l.cell2){ return cell1; }
    if(cell2!=l.cell1 && cell2!=l.cell2){ return cell2; }
    return cell3;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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()).

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 1000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1); // making a triangle loop
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else{ // strip state
       // divide one link
      CellLink dividingLink = links.get(IRand.getInt(0,links.size()-1));
      
      IVec dir = dividingLink.oppositeDir(this); // vector to the other cell
      dir.projectToPlane(nml()); // division dir is projected on normal plane
      Cell child = createChild(dir); // make a child toward the link direction
      Cell c0 = dividingLink.oppositeCell(this); // opposite cell on the link
      
      if(dividingLink.faces.size()==1){
        CellFace f1 = dividingLink.faces.get(0); // existing face on the link
        Cell c1 = f1.oppositeCell(dividingLink); // opposite cell from the link on the face
        // delete 1 link + 1 face, create 3 links + 2 faces
        new CellLink(this, child);
        new CellLink(c1,child);
        new CellLink(c0,child);
        new CellFace(this, c1, child);
        new CellFace(child, c1, c0);
        dividingLink.del();
        f1.del();
      }
      else if(dividingLink.faces.size()==2){
        CellFace f1 = dividingLink.faces.get(0);
        CellFace f2 = dividingLink.faces.get(1);
        Cell c1 = f1.oppositeCell(dividingLink);
        Cell c2 = f2.oppositeCell(dividingLink);
        // delete 1 link + 1 face, create 4 links + 4 faces
        new CellLink(this, child);
        new CellLink(c1,child);
        new CellLink(c2,child);
        new CellLink(c0,child);
        new CellFace(this, c1, child);
        new CellFace(this, child, c2);
        new CellFace(child, c1, c0);
        new CellFace(child, c0, c2);
        dividingLink.del();
        f1.del();
        f2.del();
      }
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell oppositeCell(CellLink l){ // find a cell opposite of a link edge in a triangle
    if(cell1!=l.cell1 && cell1!=l.cell2){ return cell1; }
    if(cell2!=l.cell1 && cell2!=l.cell2){ return cell2; }
    return cell3;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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.

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 1000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      // making a triangle loop
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1);
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else if(links.size()==2){ // triangle
      // making a tetrahedron enclosure
      Cell child = createChild(IRand.dir());
      CellFace f = faces.get(0);
      IVec 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
      new CellLink(f.cell1, child);
      new CellLink(f.cell2, child);
      new CellLink(f.cell3, child);
      CellFace f1 = new CellFace(f.cell1, f.cell2, child);
      CellFace f2 = new CellFace(f.cell2, f.cell3, child);
      CellFace f3 = new 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.size() > 3
      // divide one face
      CellFace dividingFace = faces.get(IRand.getInt(0,faces.size()-1));
      Cell[] c = dividingFace.oppositeCell(this);
      IVec dir = dividingFace.center().dif(pos());
      dir.projectToPlane(nml()); 
      Cell child = createChild(dir);
      // delete 1 face, create 3 links + 3 faces
      new CellLink(this, child);
      new CellLink(child, c[0]);
      new CellLink(child, c[1]);
      new CellFace(this, c[0], child);
      new CellFace(this, child, c[1]);
      new CellFace(child, c[0], c[1]);
      dividingFace.del();
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  Cell oppositeCell(CellLink l){ // find a cell opposite of a link edge in a triangle
    if(cell1!=l.cell1 && cell1!=l.cell2){ return cell1; }
    if(cell2!=l.cell1 && cell2!=l.cell2){ return cell2; }
    return cell3;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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.

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 1000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      // making a triangle loop
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1);
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else if(links.size()==2){ // triangle
      // making a tetrahedron enclosure
      Cell child = createChild(IRand.dir());
      CellFace f = faces.get(0);
      IVec 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
      new CellLink(f.cell1, child);
      new CellLink(f.cell2, child);
      new CellLink(f.cell3, child);
      CellFace f1 = new CellFace(f.cell1, f.cell2, child);
      CellFace f2 = new CellFace(f.cell2, f.cell3, child);
      CellFace f3 = new 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.size() > 3
      // divide one link
      CellLink dividingLink = links.get(IRand.getInt(0,links.size()-1));
      IVec dir = dividingLink.oppositeDir(this);
      dir.projectToPlane(nml());
      Cell child = createChild(dir);
      CellFace f1 = dividingLink.faces.get(0);
      CellFace f2 = dividingLink.faces.get(1);
      Cell c0 = dividingLink.oppositeCell(this);
      Cell c1 = f1.oppositeCell(dividingLink);
      Cell c2 = f2.oppositeCell(dividingLink);
      // delete 1 link + 2 faces, create 4 links + 4 faces
      new CellLink(this, child);
      new CellLink(c0,child);
      new CellLink(c1,child);
      new CellLink(c2,child);
      new CellFace(this, c1, child);
      new CellFace(this, child, c2);
      new CellFace(child, c1, c0);
      new CellFace(child, c0, c2);
      dividingLink.del();
      f1.del();
      f2.del();
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  Cell oppositeCell(CellLink l){ // find a cell opposite of a link edge in a triangle
    if(cell1!=l.cell1 && cell1!=l.cell2){ return cell1; }
    if(cell2!=l.cell1 && cell2!=l.cell2){ return cell2; }
    return cell3;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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.

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 1200; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
    this.push(nml().len(20)); // pressure toward normal
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      // making a triangle loop
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1);
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else if(links.size()==2){ // triangle
      // making a tetrahedron enclosure
      Cell child = createChild(IRand.dir());
      CellFace f = faces.get(0);
      IVec 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
      new CellLink(f.cell1, child);
      new CellLink(f.cell2, child);
      new CellLink(f.cell3, child);
      CellFace f1 = new CellFace(f.cell1, f.cell2, child);
      CellFace f2 = new CellFace(f.cell2, f.cell3, child);
      CellFace f3 = new 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.size() > 3
      // pick two links and insert two faces between
      int[] linkIdx = findDividingLinks();
      ArrayList< CellLink > linksToChild = new ArrayList< CellLink >(); // links to reconnect to child
      for(int i=linkIdx[0]; i <= linkIdx[1]; i++){
        linksToChild.add(links.get(i%links.size()));
      }
      IVec dir = new IVec();
      for(int i=0; i < linksToChild.size(); i++){ // average of links1 dir
        dir.add(linksToChild.get(i).oppositeDir(this).unit());
      }
      dir.projectToPlane(nml());
      Cell child = createChild(dir);
      insertChild(child, linksToChild);
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
  
  void sortLinks(){ // sort links by following connected faces
    if(links.size() <= 3) return; // no need to sort
    ArrayList< CellLink > sorted = new ArrayList< CellLink >();
    CellLink currentLink = links.get(0);
    CellFace currentFace = currentLink.faces.get(0);
    for(int i=0; i < links.size(); i++){
      sorted.add(currentLink);
      currentLink = currentFace.oppositeLink(currentLink.oppositeCell(this));
      currentFace = currentLink.oppositeFace(currentFace);
    }
    links = sorted;
  }
  
  int[] findDividingLinks(){ // find splittable two links and return the index numbers
    int linkInterval = links.size()/2; // interaval between two links 
    
    sortLinks(); // sort the order of links around cell
    int idx = IRand.getInt(0,links.size()-1);
    return new int[]{ idx, idx+linkInterval }; // index can be larger than links.size()
  }
  
  void insertChild(Cell child, ArrayList< CellLink > linksToChild){ // insert child cell and reconnect links
    ArrayList< CellFace > removingFaces = new ArrayList< CellFace >();
    for(int i=0; i < faces.size(); i++){ 
      for(int j=0; j < linksToChild.size()-1; j++){  
        if(faces.get(i).contains(linksToChild.get(j)) &&
           faces.get(i).contains(linksToChild.get(j+1))){ //a face between replacing links is removed
          removingFaces.add(faces.get(i));
          break;
        }
      }
    }
    for(int i=1; i < linksToChild.size()-1; i++){ 
      linksToChild.get(i).del(); // replacing links are once removed
      new CellLink(linksToChild.get(i).oppositeCell(this), child); // then recreated
    }
    Cell cell1 = linksToChild.get(0).oppositeCell(this); // one of two dividing link cell
    Cell cell2 = linksToChild.get(linksToChild.size()-1).oppositeCell(this); // another dividing link cell
    new CellLink(child, cell1);
    new CellLink(child, cell2);
    new CellLink(this, child);
    for(int i=0; i < removingFaces.size(); i++){
      removingFaces.get(i).del(); // replace face by deleting and recreating
      Cell[] cells = removingFaces.get(i).oppositeCell(this);
      new CellFace(child, cells[0], cells[1]);
    }
    new CellFace(this, cell1, child);
    new CellFace(this, child, cell2);
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  CellFace oppositeFace(CellFace f){ // find other face on the link
    if(faces.size()!=2) return null;
    if(faces.get(0)==f) return faces.get(1);
    return faces.get(0); 
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  CellLink oppositeLink(Cell c){ // find opposite link of a cell in a triangle
    if(!link1.contains(c)) return link1;
    if(!link2.contains(c)) return link2;
    if(!link3.contains(c)) return link3;
    return null;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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.

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 20000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  int maxLink = 6;  //limit on number of links per cell
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
    this.push(nml().len(20)); // pressure toward normal
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      // making a triangle loop
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1);
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else if(links.size()==2){ // triangle
      // making a tetrahedron enclosure
      Cell child = createChild(IRand.dir());
      CellFace f = faces.get(0);
      IVec 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
      new CellLink(f.cell1, child);
      new CellLink(f.cell2, child);
      new CellLink(f.cell3, child);
      CellFace f1 = new CellFace(f.cell1, f.cell2, child);
      CellFace f2 = new CellFace(f.cell2, f.cell3, child);
      CellFace f3 = new 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.size() > 3
      // pick two links and insert two faces between
      int[] linkIdx = findDividingLinks();
      if(linkIdx==null){ // no valid edge to split. skip division.
        active = false;
        return; 
      }
      ArrayList< CellLink > linksToChild = new ArrayList< CellLink >(); // links to reconnect to child
      for(int i=linkIdx[0]; i <= linkIdx[1]; i++){
        linksToChild.add(links.get(i%links.size()));
      }
      IVec dir = new IVec();
      for(int i=0; i < linksToChild.size(); i++){ // average of links1 dir
        dir.add(linksToChild.get(i).oppositeDir(this).unit());
      }
      dir.projectToPlane(nml());
      Cell child = createChild(dir);
      insertChild(child, linksToChild);
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
  
  void sortLinks(){ // sort links by following connected faces
    if(links.size() <= 3) return; // no need to sort
    ArrayList< CellLink > sorted = new ArrayList< CellLink >();
    CellLink currentLink = links.get(0);
    CellFace currentFace = currentLink.faces.get(0);
    for(int i=0; i < links.size(); i++){
      sorted.add(currentLink);
      currentLink = currentFace.oppositeLink(currentLink.oppositeCell(this));
      currentFace = currentLink.oppositeFace(currentFace);
    }
    links = sorted;
  }
  
  int[] findDividingLinks(){ // find splittable two links and return the index numbers
    int linkInterval = links.size()/2; // interaval between two links 
    
    sortLinks(); // sort the order of links around cell
    if(linkInterval==1 && links.size() >= maxLink){ // dividing links next to each other adds one more link
      return null;
    }
    int idx = IRand.getInt(0,links.size()-1);
    for(int i=0; i < links.size(); i++){ // check all pairs
      Cell c1 = links.get((i+idx)%links.size()).oppositeCell(this);
      Cell c2 = links.get((i+idx+linkInterval)%links.size()).oppositeCell(this);  
      if(c1.links.size() < maxLink && c2.links.size() < maxLink){ // division adds one link on the end of each dividing link
        return new int[]{ i+idx, i+idx+linkInterval }; // index can be larger than links.size()
      }
    }
    return null;
  }
  
  void insertChild(Cell child, ArrayList< CellLink > linksToChild){ // insert child cell and reconnect links
    ArrayList< CellFace > removingFaces = new ArrayList< CellFace >();
    for(int i=0; i < faces.size(); i++){ 
      for(int j=0; j < linksToChild.size()-1; j++){  
        if(faces.get(i).contains(linksToChild.get(j)) &&
           faces.get(i).contains(linksToChild.get(j+1))){ //a face between replacing links is removed
          removingFaces.add(faces.get(i));
          break;
        }
      }
    }
    for(int i=1; i < linksToChild.size()-1; i++){ 
      linksToChild.get(i).del(); // replacing links are once removed
      new CellLink(linksToChild.get(i).oppositeCell(this), child); // then recreated
    }
    Cell cell1 = linksToChild.get(0).oppositeCell(this); // one of two dividing link cell
    Cell cell2 = linksToChild.get(linksToChild.size()-1).oppositeCell(this); // another dividing link cell
    new CellLink(child, cell1);
    new CellLink(child, cell2);
    new CellLink(this, child);
    for(int i=0; i < removingFaces.size(); i++){
      removingFaces.get(i).del(); // replace face by deleting and recreating
      Cell[] cells = removingFaces.get(i).oppositeCell(this);
      new CellFace(child, cells[0], cells[1]);
    }
    new CellFace(this, cell1, child);
    new CellFace(this, child, cell2);
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  CellFace oppositeFace(CellFace f){ // find other face on the link
    if(faces.size()!=2) return null;
    if(faces.get(0)==f) return faces.get(1);
    return faces.get(0); 
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  CellLink oppositeLink(Cell c){ // find opposite link of a cell in a triangle
    if(!link1.contains(c)) return link1;
    if(!link2.contains(c)) return link2;
    if(!link3.contains(c)) return link3;
    return null;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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.

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
}

class Cell extends IParticle{
  int growthDuration = 2000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 50;
  double maxRadius = 100;
  double growthSpeed=0.1;
  int maxLink = 6; //limit on number of links per cell
  double extraLinkAllowance = 20; //percentage to accept one more link than maxLink
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
    this.push(nml().len(20)); // pressure toward normal
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
        if(IRand.pct(50)){ // random division
          active=true;
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      // making a triangle loop
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1);
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else if(links.size()==2){ // triangle
      // making a tetrahedron enclosure
      Cell child = createChild(IRand.dir());
      CellFace f = faces.get(0);
      IVec 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
      new CellLink(f.cell1, child);
      new CellLink(f.cell2, child);
      new CellLink(f.cell3, child);
      CellFace f1 = new CellFace(f.cell1, f.cell2, child);
      CellFace f2 = new CellFace(f.cell2, f.cell3, child);
      CellFace f3 = new 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.size() > 3
      // pick two links and insert two faces between
      int[] linkIdx = findDividingLinks();
      if(linkIdx==null){ // no valid edge to split. skip division.
        active = false;
        return; 
      }
      ArrayList< CellLink > linksToChild = new ArrayList< CellLink >(); // links to reconnect to child
      for(int i=linkIdx[0]; i <= linkIdx[1]; i++){
        linksToChild.add(links.get(i%links.size()));
      }
      IVec dir = new IVec();
      for(int i=0; i < linksToChild.size(); i++){ // average of links1 dir
        dir.add(linksToChild.get(i).oppositeDir(this).unit());
      }
      dir.projectToPlane(nml());
      Cell child = createChild(dir);
      insertChild(child, linksToChild);
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false;  //reset activation
    return child;
  }
  
  void sortLinks(){ // sort links by following connected faces
    if(links.size() <= 3) return; // no need to sort
    ArrayList< CellLink > sorted = new ArrayList< CellLink >();
    CellLink currentLink = links.get(0);
    CellFace currentFace = currentLink.faces.get(0);
    for(int i=0; i < links.size(); i++){
      sorted.add(currentLink);
      currentLink = currentFace.oppositeLink(currentLink.oppositeCell(this));
      currentFace = currentLink.oppositeFace(currentFace);
    }
    links = sorted;
  }
  
  int[] findDividingLinks(){ // find splittable two links and return the index numbers
    int linkInterval = links.size()/2; // interaval between two links 
    
    sortLinks(); // sort the order of links around cell
    boolean allowExtra = IRand.pct(extraLinkAllowance);
    if(linkInterval==1 && // dividing links next to each other adds one more link
       !(links.size() < maxLink || allowExtra && links.size() < maxLink+1)){
      return null;
    }
    int idx = IRand.getInt(0,links.size()-1);
    for(int i=0; i < links.size(); i++){ // check all pairs
      Cell c1 = links.get((i+idx)%links.size()).oppositeCell(this);
      Cell c2 = links.get((i+idx+linkInterval)%links.size()).oppositeCell(this);  
      if(c1.links.size() < maxLink && c2.links.size() < maxLink || // division adds one link on the end of each dividing link
         allowExtra && c1.links.size() < maxLink+1 && c2.links.size() < maxLink+1){
        return new int[]{ i+idx, i+idx+linkInterval }; // index can be larger than links.size()
      }
    }
    return null;
  }
  
  void insertChild(Cell child, ArrayList< CellLink > linksToChild){ // insert child cell and reconnect links
    ArrayList< CellFace > removingFaces = new ArrayList< CellFace >();
    for(int i=0; i < faces.size(); i++){ 
      for(int j=0; j < linksToChild.size()-1; j++){  
        if(faces.get(i).contains(linksToChild.get(j)) &&
           faces.get(i).contains(linksToChild.get(j+1))){ //a face between replacing links is removed
          removingFaces.add(faces.get(i));
          break;
        }
      }
    }
    for(int i=1; i < linksToChild.size()-1; i++){ 
      linksToChild.get(i).del(); // replacing links are once removed
      new CellLink(linksToChild.get(i).oppositeCell(this), child); // then recreated
    }
    Cell cell1 = linksToChild.get(0).oppositeCell(this); // one of two dividing link cell
    Cell cell2 = linksToChild.get(linksToChild.size()-1).oppositeCell(this); // another dividing link cell
    new CellLink(child, cell1);
    new CellLink(child, cell2);
    new CellLink(this, child);
    for(int i=0; i < removingFaces.size(); i++){
      removingFaces.get(i).del(); // replace face by deleting and recreating
      Cell[] cells = removingFaces.get(i).oppositeCell(this);
      new CellFace(child, cells[0], cells[1]);
    }
    new CellFace(this, cell1, child);
    new CellFace(this, child, cell2);
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  CellFace oppositeFace(CellFace f){ // find other face on the link
    if(faces.size()!=2) return null;
    if(faces.get(0)==f) return faces.get(1);
    return faces.get(0); 
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  CellLink oppositeLink(Cell c){ // find opposite link of a cell in a triangle
    if(!link1.contains(c)) return link1;
    if(!link2.contains(c)) return link2;
    if(!link3.contains(c)) return link3;
    return null;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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.

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
  cell.active = true;
}

class Cell extends IParticle{
  int growthDuration = 10000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 20;
  double maxRadius = 100;
  double growthSpeed=0.1;
  int maxLink = 6; //limit on number of links per cell
  double extraLinkAllowance = 0; //percentage to accept one more link than maxLink
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(100); // constant force
      this.push(dif);
    }
    this.push(nml().len(20)); // pressure toward normal
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      // making a triangle loop
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1);
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else if(links.size()==2){ // triangle
      // making a tetrahedron enclosure
      Cell child = createChild(IRand.dir());
      CellFace f = faces.get(0);
      IVec 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
      new CellLink(f.cell1, child);
      new CellLink(f.cell2, child);
      new CellLink(f.cell3, child);
      CellFace f1 = new CellFace(f.cell1, f.cell2, child);
      CellFace f2 = new CellFace(f.cell2, f.cell3, child);
      CellFace f3 = new 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.size() > 3
      // pick two links and insert two faces between
      int[] linkIdx = findDividingLinks();
      if(linkIdx==null){ // no valid edge to split. skip division.
        active = false;
        return; 
      }
      ArrayList< CellLink > linksToChild = new ArrayList< CellLink >(); // links to reconnect to child
      for(int i=linkIdx[0]; i <= linkIdx[1]; i++){
        linksToChild.add(links.get(i%links.size()));
      }
      IVec dir = new IVec();
      for(int i=0; i < linksToChild.size(); i++){ // average of links1 dir
        dir.add(linksToChild.get(i).oppositeDir(this).unit());
      }
      dir.projectToPlane(nml());
      Cell child = createChild(dir);
      insertChild(child, linksToChild);
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    active=false; //reset activation
    child.active = true; //activate child always
    return child;
  }
  
  void sortLinks(){ // sort links by following connected faces
    if(links.size() <= 3) return; // no need to sort
    ArrayList< CellLink > sorted = new ArrayList< CellLink >();
    CellLink currentLink = links.get(0);
    CellFace currentFace = currentLink.faces.get(0);
    for(int i=0; i < links.size(); i++){
      sorted.add(currentLink);
      currentLink = currentFace.oppositeLink(currentLink.oppositeCell(this));
      currentFace = currentLink.oppositeFace(currentFace);
    }
    links = sorted;
  }
  
  int[] findDividingLinks(){ // find splittable two links and return the index numbers
    int linkInterval = 1; // face division only 
    
    sortLinks(); // sort the order of links around cell
    boolean allowExtra = IRand.pct(extraLinkAllowance);
    if(linkInterval==1 && // dividing links next to each other adds one more link
       !(links.size() < maxLink || allowExtra && links.size() < maxLink+1)){
      return null;
    }
    int idx = IRand.getInt(0,links.size()-1);
    for(int i=0; i < links.size(); i++){ // check all pairs
      Cell c1 = links.get((i+idx)%links.size()).oppositeCell(this);
      Cell c2 = links.get((i+idx+linkInterval)%links.size()).oppositeCell(this);  
      if(c1.links.size() < maxLink && c2.links.size() < maxLink || // division adds one link on the end of each dividing link
         allowExtra && c1.links.size() < maxLink+1 && c2.links.size() < maxLink+1){
        return new int[]{ i+idx, i+idx+linkInterval }; // index can be larger than links.size()
      }
    }
    return null;
  }
  
  void insertChild(Cell child, ArrayList< CellLink > linksToChild){ // insert child cell and reconnect links
    ArrayList< CellFace > removingFaces = new ArrayList< CellFace >();
    for(int i=0; i < faces.size(); i++){ 
      for(int j=0; j < linksToChild.size()-1; j++){  
        if(faces.get(i).contains(linksToChild.get(j)) &&
           faces.get(i).contains(linksToChild.get(j+1))){ //a face between replacing links is removed
          removingFaces.add(faces.get(i));
          break;
        }
      }
    }
    for(int i=1; i < linksToChild.size()-1; i++){ 
      linksToChild.get(i).del(); // replacing links are once removed
      new CellLink(linksToChild.get(i).oppositeCell(this), child); // then recreated
    }
    Cell cell1 = linksToChild.get(0).oppositeCell(this); // one of two dividing link cell
    Cell cell2 = linksToChild.get(linksToChild.size()-1).oppositeCell(this); // another dividing link cell
    new CellLink(child, cell1);
    new CellLink(child, cell2);
    new CellLink(this, child);
    for(int i=0; i < removingFaces.size(); i++){
      removingFaces.get(i).del(); // replace face by deleting and recreating
      Cell[] cells = removingFaces.get(i).oppositeCell(this);
      new CellFace(child, cells[0], cells[1]);
    }
    new CellFace(this, cell1, child);
    new CellFace(this, child, cell2);
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  CellFace oppositeFace(CellFace f){ // find other face on the link
    if(faces.size()!=2) return null;
    if(faces.get(0)==f) return faces.get(1);
    return faces.get(0); 
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  CellLink oppositeLink(Cell c){ // find opposite link of a cell in a triangle
    if(!link1.contains(c)) return link1;
    if(!link2.contains(c)) return link2;
    if(!link3.contains(c)) return link3;
    return null;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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.

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
  cell.active = true;
}

class Cell extends IParticle{
  int growthDuration = 2000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 20;
  double maxRadius = 100;
  double growthSpeed=0.1;
  int maxLink = 8; //limit on number of links per cell
  double extraLinkAllowance = 0; //percentage to accept one more link than maxLink
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(100); // constant force
      this.push(dif);
    }
    this.push(nml().len(20)); // pressure toward normal
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      // making a triangle loop
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1);
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else if(links.size()==2){ // triangle
      // making a tetrahedron enclosure
      Cell child = createChild(IRand.dir());
      CellFace f = faces.get(0);
      IVec 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
      new CellLink(f.cell1, child);
      new CellLink(f.cell2, child);
      new CellLink(f.cell3, child);
      CellFace f1 = new CellFace(f.cell1, f.cell2, child);
      CellFace f2 = new CellFace(f.cell2, f.cell3, child);
      CellFace f3 = new 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.size() > 3
      // pick two links and insert two faces between
      int[] linkIdx = findDividingLinks();
      if(linkIdx==null){ // no valid edge to split. skip division.
        active = false;
        return; 
      }
      ArrayList< CellLink > linksToChild = new ArrayList< CellLink >(); // links to reconnect to child
      for(int i=linkIdx[0]; i <= linkIdx[1]; i++){
        linksToChild.add(links.get(i%links.size()));
      }
      IVec dir = new IVec();
      for(int i=0; i < linksToChild.size(); i++){ // average of links1 dir
        dir.add(linksToChild.get(i).oppositeDir(this).unit());
      }
      dir.projectToPlane(nml());
      Cell child = createChild(dir);
      insertChild(child, linksToChild);
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    if(IRand.pct(90)){ active=false; }//10% stay active
    child.active = true; //activate child always
    return child;
  }
  
  void sortLinks(){ // sort links by following connected faces
    if(links.size() <= 3) return; // no need to sort
    ArrayList< CellLink > sorted = new ArrayList< CellLink >();
    CellLink currentLink = links.get(0);
    CellFace currentFace = currentLink.faces.get(0);
    for(int i=0; i < links.size(); i++){
      sorted.add(currentLink);
      currentLink = currentFace.oppositeLink(currentLink.oppositeCell(this));
      currentFace = currentLink.oppositeFace(currentFace);
    }
    links = sorted;
  }
  
  int[] findDividingLinks(){ // find splittable two links and return the index numbers
    int linkInterval = 1; // interaval between two links 
    
    sortLinks(); // sort the order of links around cell
    boolean allowExtra = IRand.pct(extraLinkAllowance);
    if(linkInterval==1 && // dividing links next to each other adds one more link
       !(links.size() < maxLink || allowExtra && links.size() < maxLink+1)){
      return null;
    }
    int idx = IRand.getInt(0,links.size()-1);
    for(int i=0; i < links.size(); i++){ // check all pairs
      Cell c1 = links.get((i+idx)%links.size()).oppositeCell(this);
      Cell c2 = links.get((i+idx+linkInterval)%links.size()).oppositeCell(this);  
      if(c1.links.size() < maxLink && c2.links.size() < maxLink || // division adds one link on the end of each dividing link
         allowExtra && c1.links.size() < maxLink+1 && c2.links.size() < maxLink+1){
        return new int[]{ i+idx, i+idx+linkInterval }; // index can be larger than links.size()
      }
    }
    return null;
  }
  
  void insertChild(Cell child, ArrayList< CellLink > linksToChild){ // insert child cell and reconnect links
    ArrayList< CellFace > removingFaces = new ArrayList< CellFace >();
    for(int i=0; i < faces.size(); i++){ 
      for(int j=0; j < linksToChild.size()-1; j++){  
        if(faces.get(i).contains(linksToChild.get(j)) &&
           faces.get(i).contains(linksToChild.get(j+1))){ //a face between replacing links is removed
          removingFaces.add(faces.get(i));
          break;
        }
      }
    }
    for(int i=1; i < linksToChild.size()-1; i++){ 
      linksToChild.get(i).del(); // replacing links are once removed
      new CellLink(linksToChild.get(i).oppositeCell(this), child); // then recreated
    }
    Cell cell1 = linksToChild.get(0).oppositeCell(this); // one of two dividing link cell
    Cell cell2 = linksToChild.get(linksToChild.size()-1).oppositeCell(this); // another dividing link cell
    new CellLink(child, cell1);
    new CellLink(child, cell2);
    new CellLink(this, child);
    for(int i=0; i < removingFaces.size(); i++){
      removingFaces.get(i).del(); // replace face by deleting and recreating
      Cell[] cells = removingFaces.get(i).oppositeCell(this);
      new CellFace(child, cells[0], cells[1]);
    }
    new CellFace(this, cell1, child);
    new CellFace(this, child, cell2);
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  CellFace oppositeFace(CellFace f){ // find other face on the link
    if(faces.size()!=2) return null;
    if(faces.get(0)==f) return faces.get(1);
    return faces.get(0); 
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  CellLink oppositeLink(Cell c){ // find opposite link of a cell in a triangle
    if(!link1.contains(c)) return link1;
    if(!link2.contains(c)) return link2;
    if(!link3.contains(c)) return link3;
    return null;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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.

import igeo.*;
import processing.opengl.*;

IMesh mesh; // global variable of mesh geometry

void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  mesh = new IMesh().clr(0.7,0,0);
  Cell cell = new Cell(new IVec(0,0,0), 10);
  cell.stayActive = true;
  cell.activateChild = true;
  cell.active = true;
}

class Cell extends IParticle{
  int growthDuration = 1800; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 20;
  double maxRadius = 100;
  double growthSpeed=0.1;
  int maxLink = 6; //limit on number of links per cell
  double extraLinkAllowance = 10; //percentage to accept one more link than maxLink
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  boolean stayActive = false; //activation control variables
  boolean activateChild = false;
  IVertex vertex;
  
  Cell(IVec pos, double rad){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len(((c.radius+radius)-dif.len())*100+50); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter).len(50); // constant force
      this.push(dif);
    }
    this.push(nml().len(100)); // pressure toward normal
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()==0){ // dot state
      Cell child = createChild(IRand.dir());
      new CellLink(this, child); // add one link
    }
    else if(links.size()==1){ // line state 
      // making a triangle loop
      Cell child = createChild(IRand.dir());
      new CellLink(child, links.get(0).cell1);
      new CellLink(child, links.get(0).cell2);
      new CellFace(child, links.get(0).cell1, links.get(0).cell2); // making a triangle face
    }
    else if(links.size()==2){ // triangle
      // making a tetrahedron enclosure
      Cell child = createChild(IRand.dir());
      CellFace f = faces.get(0);
      IVec 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
      new CellLink(f.cell1, child);
      new CellLink(f.cell2, child);
      new CellLink(f.cell3, child);
      CellFace f1 = new CellFace(f.cell1, f.cell2, child);
      CellFace f2 = new CellFace(f.cell2, f.cell3, child);
      CellFace f3 = new 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.size() > 3
      // pick two links and insert two faces between
      int[] linkIdx = findDividingLinks();
      if(linkIdx==null){ // no valid edge to split. skip division.
        active = false;
        return; 
      }
      ArrayList< CellLink > linksToChild = new ArrayList< CellLink >(); // links to reconnect to child
      for(int i=linkIdx[0]; i <= linkIdx[1]; i++){
        linksToChild.add(links.get(i%links.size()));
      }
      IVec dir = new IVec();
      for(int i=0; i < linksToChild.size(); i++){ // average of links1 dir
        dir.add(linksToChild.get(i).oppositeDir(this).unit());
      }
      dir.projectToPlane(nml());
      Cell child = createChild(dir);
      insertChild(child, linksToChild);
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius);
    pos().sub(dir); // move to the other way
    if(stayActive){ active=true; } //activation control
    else{ active=false; } 
    if(activateChild){ child.active=true; } 
    else{ child.active=false; }
    child.stayActive = stayActive; //activation behavior parameters are passed to child
    child.activateChild = activateChild;
    return child;
  }
  
  void sortLinks(){ // sort links by following connected faces
    if(links.size() <= 3) return; // no need to sort
    ArrayList< CellLink > sorted = new ArrayList< CellLink >();
    CellLink currentLink = links.get(0);
    CellFace currentFace = currentLink.faces.get(0);
    for(int i=0; i < links.size(); i++){
      sorted.add(currentLink);
      currentLink = currentFace.oppositeLink(currentLink.oppositeCell(this));
      currentFace = currentLink.oppositeFace(currentFace);
    }
    links = sorted;
  }
  
  int[] findDividingLinks(){ // find splittable two links and return the index numbers
    int linkInterval = links.size()/2; // interaval between two links 
    
    sortLinks(); // sort the order of links around cell
    boolean allowExtra = IRand.pct(extraLinkAllowance);
    if(linkInterval==1 && // dividing links next to each other adds one more link
       !(links.size() < maxLink || allowExtra && links.size() < maxLink+1)){
      return null;
    }
    int idx = IRand.getInt(0,links.size()-1);
    for(int i=0; i < links.size(); i++){ // check all pairs
      Cell c1 = links.get((i+idx)%links.size()).oppositeCell(this);
      Cell c2 = links.get((i+idx+linkInterval)%links.size()).oppositeCell(this);  
      if(c1.links.size() < maxLink && c2.links.size() < maxLink || // division adds one link on the end of each dividing link
         allowExtra && c1.links.size() < maxLink+1 && c2.links.size() < maxLink+1){
        return new int[]{ i+idx, i+idx+linkInterval }; // index can be larger than links.size()
      }
    }
    return null;
  }
  
  void insertChild(Cell child, ArrayList< CellLink > linksToChild){ // insert child cell and reconnect links
    ArrayList< CellFace > removingFaces = new ArrayList< CellFace >();
    for(int i=0; i < faces.size(); i++){ 
      for(int j=0; j < linksToChild.size()-1; j++){  
        if(faces.get(i).contains(linksToChild.get(j)) &&
           faces.get(i).contains(linksToChild.get(j+1))){ //a face between replacing links is removed
          removingFaces.add(faces.get(i));
          break;
        }
      }
    }
    for(int i=1; i < linksToChild.size()-1; i++){ 
      linksToChild.get(i).del(); // replacing links are once removed
      new CellLink(linksToChild.get(i).oppositeCell(this), child); // then recreated
    }
    Cell cell1 = linksToChild.get(0).oppositeCell(this); // one of two dividing link cell
    Cell cell2 = linksToChild.get(linksToChild.size()-1).oppositeCell(this); // another dividing link cell
    new CellLink(child, cell1);
    new CellLink(child, cell2);
    new CellLink(this, child);
    for(int i=0; i < removingFaces.size(); i++){
      removingFaces.get(i).del(); // replace face by deleting and recreating
      Cell[] cells = removingFaces.get(i).oppositeCell(this);
      new CellFace(child, cells[0], cells[1]);
    }
    new CellFace(this, cell1, child);
    new CellFace(this, child, cell2);
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  CellFace oppositeFace(CellFace f){ // find other face on the link
    if(faces.size()!=2) return null;
    if(faces.get(0)==f) return faces.get(1);
    return faces.get(0); 
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  
  CellFace(Cell c1, Cell c2, Cell c3){
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  CellLink oppositeLink(Cell c){ // find opposite link of a cell in a triangle
    if(!link1.contains(c)) return link1;
    if(!link2.contains(c)) return link2;
    if(!link3.contains(c)) return link3;
    return null;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(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

import igeo.*;
import processing.opengl.*;

double averageLength; // to adjust cell radius based on size
double sizeRatio; // to adjustment forces based on size
 
void setup(){
  size(480,360,IG.GL);
  IConfig.syncDrawAndDynamics=true; //not to crash when some geometry is deleted while drawing
  IG.open("mesh2.3dm");
  
  IPoint[] fixedPts = IG.layer("fix points").points(); // points to speficy fixed cells
  IPoint[] activePts = IG.layer("active points").points(); // points to speficy active cells
  for(int i=0; i < IG.meshNum(); i++){
    IMesh m = IG.mesh(i);
    m.triangulate(); // quad mesh not accepted
    
    Cell[] cells = cellsFromMesh(m, fixedPts);
    
    for(int j=0; j < cells.length; j++){
      for(int k=0; activePts!=null && k < activePts.length; k++){
        if(cells[j].pos().eq(activePts[k])){ // activate cell if active point is located
          cells[j].stayActiveProbability = 80;
          cells[j].activateChild=true;
          cells[j].active=true;
        }
      }
    }
  }
  IG.del(activePts);
}

Cell[] cellsFromMesh(IMesh mesh, IPoint[] fixedPoints){ // create cells, cell links, cell faces from a mesh
  Cell[] cells = new Cell[mesh.vertexNum()];
  averageLength=0;
  for(int i=0; i < mesh.edgeNum(); i++){
    averageLength += mesh.edge(i).vertex(0).dist(mesh.edge(i).vertex(1));
  }
  averageLength /= mesh.edgeNum();
  sizeRatio = averageLength/10;
  
  for(int i=0; i < mesh.vertexNum(); i++){
    double radius = 0;
    for(int j=0; j < mesh.vertex(i).edgeNum(); j++){
      radius += mesh.vertex(i).edge(j).vertex(0).dist(mesh.vertex(i).edge(j).vertex(1)) / 2;
    }
    radius /= mesh.vertex(i).edgeNum(); // average
    cells[i] = new Cell(mesh.vertex(i), radius, mesh);
    
    boolean done=false; 
    for(int j=0; j < fixedPoints.length&&!done; j++){
      if(cells[i].pos.eq(fixedPoints[j])){
        cells[i].fix();
      }
    }
  }
  
  CellLink[] links = new CellLink[mesh.edgeNum()];
  for(int i=0; i < mesh.edgeNum(); i++){
    int idx1 = mesh.mesh.vertices.indexOf(mesh.edge(i).vertex(0));
    int idx2 = mesh.mesh.vertices.indexOf(mesh.edge(i).vertex(1));
    links[i] = new CellLink(cells[idx1], cells[idx2]);
  }
  
  CellFace[] faces = new CellFace[mesh.faceNum()];
  for(int i=0; i < mesh.faceNum(); i++){
    int idx1 = mesh.mesh.vertices.indexOf(mesh.face(i).vertex(0));
    int idx2 = mesh.mesh.vertices.indexOf(mesh.face(i).vertex(1));
    int idx3 = mesh.mesh.vertices.indexOf(mesh.face(i).vertex(2));
    faces[i] = new CellFace(cells[idx1],cells[idx2],cells[idx3],mesh.face(i),mesh);
  }
  
  return cells;
}

class Cell extends IParticle{
  int growthDuration = 2000; //duration of growth and division
  int growthInterval = 10;
  int divisionInterval = 100;
  double maxRadius = averageLength/2; // no bigger than half of average edge length
  double growthSpeed=0.5*sizeRatio;
  int maxLink = 6; //limit on number of links per cell
  double extraLinkAllowance = 100; //percentage to accept one more link than maxLink
  
  ArrayList< CellLink > links; //store links
  ArrayList< CellFace > faces; //store faces
  double radius;
  boolean active = false;
  double stayActiveProbability = 0; //activation control probability
  boolean activateChild = false;
  IVertex vertex;
  
  IMesh mesh;
  
  Cell(IVec pos, double rad, IMesh m){
    super(pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = new IVertex();
    vertex.pos = pos; // set reference of particle pos directly
    mesh = m;
    fric(0.2);
  }
  
  Cell(IVertex vtx, double rad, IMesh m){
    super(vtx.pos, new IVec(0,0,0));
    radius = rad;
    links = new ArrayList< CellLink >();
    faces = new ArrayList< CellFace >();
    vertex = vtx;
    mesh = m;
    fric(0.2);
  }
  
  void interact(ArrayList< IDynamics > agents){
    IVec neighborCenter = new IVec(0,0,0);
    int neighborCount=0;
    double neighborDist = radius*4; 
    for(int i=0; i < agents.size(); i++){
      if(agents.get(i) instanceof Cell){
        Cell c = (Cell)agents.get(i);
        if(c != this){
          // push if closer than two radii
          if(c.pos().dist(pos()) < c.radius+radius){
            IVec dif = c.pos().dif(pos());
            dif.len((((c.radius+radius)-dif.len())*100+50)*sizeRatio); // the closer the harder to push
            c.push(dif);
          }
          // count neighbors and calculate their center
          if(c.pos().dist(pos()) < c.radius+radius + neighborDist){
            neighborCenter.add(c.pos());
            neighborCount++;
          }
        }
      }
    }
    if(neighborCount >= 1){ // push from center of neighbors
      neighborCenter.div(neighborCount);
      IVec dif = pos().dif(neighborCenter);
      if(dif.len()>0){
        this.push(dif.len(50*sizeRatio)); // constant force
      }
    }
    this.push(nml().len(50*sizeRatio)); // pressure toward normal
  }
  
  void update(){
    if(IG.time() < growthDuration){
      if(time() > 0 && time()%divisionInterval==0){
        if(active){ // divide when active flag is on
          divide();
        }
      }
      if(time()%growthInterval==0){
        grow();
      }
      if(active){ clr(1.0,0.5,0); } // update color by being active
      else{ clr(mesh.clr()); }
    }
    vertex.nml(nml()); // update mesh vertex normal
  }
  
  void grow(){ // growing cell size
    if(radius < maxRadius){
      radius += growthSpeed;
    }  
  }
  
  void divide(){ // cell division
    if(links.size()>=3){  // divide only if it has more than 3 links
      // pick two links and insert two faces between
      int[] linkIdx = findDividingLinks();
      if(linkIdx==null){ // no valid edge to split. skip division.
        active = false;
        return; 
      }
      ArrayList< CellLink > linksToChild = new ArrayList< CellLink >(); // links to reconnect to child
      
      if(linkIdx.length==1){ // dividing open edge
        for(int i=linkIdx[0]; i < links.size(); i++){
          linksToChild.add(links.get(i%links.size()));
        }
      }
      else{
        for(int i=linkIdx[0]; i <= linkIdx[1]; i++){
          linksToChild.add(links.get(i%links.size()));
        }
      }
      IVec dir = new IVec();
      for(int i=0; i < linksToChild.size(); i++){ // average of links1 dir
        dir.add(linksToChild.get(i).oppositeDir(this).unit());
      }
      dir.projectToPlane(nml());
      Cell child = createChild(dir);
      insertChild(child, linksToChild, linkIdx.length==1);
    }
  }
  
  IVec nml(){ // calc vertex normal from face normal
    if(faces.size()==0){ return new IVec(0,0,1); }
    IVec n = new IVec(0,0,0);
    for(int i=0; i < faces.size(); i++){ n.add(faces.get(i).nml()); }
    return n.unit();
  }
  
  Cell createChild(IVec dir){
    radius *= 0.5; //make both cell size half
    dir.len(radius);
    Cell child = new Cell(pos().cp(dir), radius, mesh);
    pos().sub(dir); // move to the other way
    if(IRand.pct(stayActiveProbability)){ active=true; } //activation control
    else{ active=false; } 
    if(activateChild){ child.active=true; } 
    else{ child.active=false; }
    child.stayActiveProbability = stayActiveProbability; //activation behavior parameters are passed to child
    child.activateChild = activateChild;
    return child;
  }
  
  void sortLinks(){ // sort links by following connected faces
    if(links.size() <= 2) return; // no need to sort
    ArrayList< CellLink > sorted = new ArrayList< CellLink >();
    CellLink currentLink = null;
    for(int i=0; i < links.size() && currentLink==null; i++){
      if(links.get(i).faces.size()==1){ // open link
        currentLink = links.get(i);
      }
    }
    if(currentLink == null){ // no open link
      if(links.size()==3) return;
      currentLink = links.get(0);
    }
    
    CellFace currentFace = currentLink.faces.get(0);
    for(int i=0; i < links.size(); i++){
      sorted.add(currentLink);
      if(i < links.size()-1){
        if(currentFace==null) break;
        currentLink = currentFace.oppositeLink(currentLink.oppositeCell(this));
        if(currentLink==null) break;
        currentFace = currentLink.oppositeFace(currentFace);
      }
    }
    links = sorted;
  }
  
  int[] findDividingLinks(){ // find splittable two links and return the index numbers
    int linkInterval = links.size()/2; // interaval between two links
    
    sortLinks(); // sort the order of links around cell
    boolean allowExtra = IRand.pct(extraLinkAllowance);
    if(linkInterval==1 && // dividing links next to each other adds one more link
       !(links.size() < maxLink || allowExtra && links.size() < maxLink+1)){
      return null;
    }
    
    if(links.get(0).faces.size()==1 && IRand.pct(10)){ // if open edge, divide at one edge
      int idx = IRand.getInt(0,links.size()-1);
      for(int i=0; i < links.size(); i++){ // check all pairs
        if( ((i+idx)%links.size())!=0 && ((i+idx)%links.size())!=links.size()-1 ){ 
          Cell c1 = links.get((i+idx)%links.size()).oppositeCell(this);
          if(c1.links.size() < maxLink  || // division adds one link on the end of each dividing link
            allowExtra && c1.links.size() < maxLink+1){
            return new int[]{ (i+idx)%links.size() }; // index can be larger than links.size()
          }
        }
      }
      return null;
    }
    
    int idx = IRand.getInt(0,links.size()-1);
    for(int i=0; i < links.size(); i++){ // check all pairs
      Cell c1 = links.get((i+idx)%links.size()).oppositeCell(this);
      Cell c2 = links.get((i+idx+linkInterval)%links.size()).oppositeCell(this);  
      if(c1.links.size() < maxLink && c2.links.size() < maxLink || // division adds one link on the end of each dividing link
         allowExtra && c1.links.size() < maxLink+1 && c2.links.size() < maxLink+1){
        return new int[]{ i+idx, i+idx+linkInterval }; // index can be larger than links.size()
      }
    }
    return null;
  }
  
  void insertChild(Cell child, ArrayList< CellLink > linksToChild, boolean insertAtEdge){ // insert child cell and reconnect links
    ArrayList< CellFace > removingFaces = new ArrayList< CellFace >();
    for(int j=0; j < linksToChild.size()-1; j++){
      for(int i=0; i < faces.size(); i++){ 
        if(faces.get(i).contains(linksToChild.get(j)) &&
           faces.get(i).contains(linksToChild.get(j+1))){ //a face between replacing links is removed
          removingFaces.add(faces.get(i));
          break;
        }
      }
    }
    
    if(insertAtEdge){
      for(int i=1; i < linksToChild.size(); i++){ 
        linksToChild.get(i).del(); // replacing links are once removed
        new CellLink(linksToChild.get(i).oppositeCell(this), child); // then recreated
      }
      Cell cell1 = linksToChild.get(0).oppositeCell(this); // one of two dividing link cell
      new CellLink(child, cell1);
      new CellLink(this, child);
      for(int i=0; i < removingFaces.size(); i++){
        removingFaces.get(i).del(); // replace face by deleting and recreating
        Cell[] cells = removingFaces.get(i).oppositeCell(this);
        new CellFace(child, cells[0], cells[1], mesh);
      }
      new CellFace(this, cell1, child, mesh);
    }
    else{
      for(int i=1; i < linksToChild.size()-1; i++){ 
        linksToChild.get(i).del(); // replacing links are once removed
        new CellLink(linksToChild.get(i).oppositeCell(this), child); // then recreated
      }
      Cell cell1 = linksToChild.get(0).oppositeCell(this); // one of two dividing link cell
      Cell cell2 = linksToChild.get(linksToChild.size()-1).oppositeCell(this); // another dividing link cell
      new CellLink(child, cell1);
      new CellLink(child, cell2);
      new CellLink(this, child);
      
      for(int i=0; i < removingFaces.size(); i++){
        removingFaces.get(i).del(); // replace face by deleting and recreating
        Cell[] cells = removingFaces.get(i).oppositeCell(this);
        new CellFace(child, cells[0], cells[1], mesh);
      }
      new CellFace(this, cell1, child, mesh);
      new CellFace(this, child, cell2, mesh);
    }
  }
  
  CellLink linkTo(Cell c){ // find a link between 2 cells
    for(int i=0; i < links.size() ; i++){
      if(links.get(i).contains(c)){ return links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
}

class CellLink extends IAgent{ // a link to connect 2 cells with spring force
  double maxForce = 100;
  Cell cell1, cell2;
  ArrayList< CellFace > faces;
  
  CellLink(Cell c1, Cell c2){
    cell1 = c1; cell2 = c2;
    cell1.links.add(this); // register this link to cells{
    cell2.links.add(this);
    faces = new ArrayList< CellFace >();
  }
  
  void interact(ArrayList< IDynamics > agents){
    // spring force
    IVec dif = cell1.pos().dif(cell2.pos());
    double force = (dif.len()-(cell1.radius+cell2.radius))/(cell1.radius+cell2.radius)*300*sizeRatio;
    if(force > maxForce){ force=maxForce; } // avoid putting too much force
    else if(force < -maxForce){ force=-maxForce; }
    dif.len(force);
    cell1.pull(dif);
    cell2.push(dif);
  }
  
  boolean contains(Cell c){ //check if link contains the cell
    if(c==cell1 || c==cell2) return true;
    return false;
  }
  
  void del(){
    cell1.links.remove(this); // unregister from cells
    cell2.links.remove(this);
    super.del(); // stop agent
  }
  
  CellFace oppositeFace(CellFace f){ // find other face on the link
    if(faces.size()!=2) return null;
    if(faces.get(0)==f) return faces.get(1);
    return faces.get(0); 
  }
  
  Cell oppositeCell(Cell c){ // find other cell on the link
    if(cell1==c) return cell2;
    if(cell2==c) return cell1;
    IG.err("Link does not contain the input cell");
    return null; 
  }
  
  IVec oppositeDir(Cell c){ //calculate a vector to the other cell
    return oppositeCell(c).pos().dif(c.pos());
  }
  
  Cell sharedCell(CellLink link){
    if(link.cell1==cell1) return cell1;
    if(link.cell2==cell1) return cell1;
    if(link.cell1==cell2) return cell2;
    if(link.cell2==cell2) return cell2;
    return null; 
  }
}

class CellFace{ // a triangle grouping 3 cells and 3 links
  Cell cell1, cell2, cell3;
  CellLink link1, link2, link3;
  IFace face;
  IMesh mesh;
  
  CellFace(Cell c1, Cell c2, Cell c3, IMesh m){
    mesh = m;
    link1 = findLink(c1, c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = findLink(cell1, cell2);
    link2 = findLink(cell2, cell3);
    link3 = findLink(cell3, cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }
  
  CellFace(CellLink l1, CellLink l2, CellLink l3, IMesh m){
    mesh = m;
    Cell c1 = l3.sharedCell(l1);
    Cell c2 = l1.sharedCell(l2);
    Cell c3 = l2.sharedCell(l3);
    
    // keep order of cells in right hand screw for consistent normal
    if(l1.faces.size()==1 && l1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
      link1 = l1; link2 = l3; link3 = l2;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
      link1 = l1; link2 = l2; link3 = l3;
    }
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = new IFace(cell1.vertex,cell2.vertex,cell3.vertex);
    mesh.addFace(face); // add triangular face to the mesh
  }

  CellFace(Cell c1, Cell c2, Cell c3, IFace f, IMesh m){
    mesh = m;
    link1 = c1.linkTo(c2);
    // keep order of cells in right hand screw for consistent normal
    if(link1.faces.size()==1 && link1.faces.get(0).cellOrder(c1,c2)){
      cell1 = c2; cell2 = c1; cell3 = c3;
    }
    else{
      cell1 = c1; cell2 = c2; cell3 = c3;
    }
    link1 = cell1.linkTo(cell2);
    link2 = cell2.linkTo(cell3);
    link3 = cell3.linkTo(cell1);
    cell1.faces.add(this); // register this face to cells and links
    cell2.faces.add(this);
    cell3.faces.add(this);
    link1.faces.add(this);
    link2.faces.add(this);
    link3.faces.add(this);
    face = f;
  }

  
  IVec center(){ // calc center of triangle
    return IVec.center(cell1.pos(), cell2.pos(), cell3.pos());
  }
  
  IVec nml(){ // normal vector
    return cell1.pos().nml(cell2.pos(), cell3.pos()).unit();
  }
  
  void flipNml(){ // flip normal
    Cell tmp = cell2;
    CellLink tmpLink = link2;
    cell2 = cell3;
    link2 = link3;
    cell3 = tmp;
    link3 = tmpLink;
  }
  
  // true if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise false
  boolean cellOrder(Cell c1, Cell c2){
    if(cell1==c1 && cell2==c2 || cell2==c1 && cell3==c2 || cell3==c1 && cell1==c2 ){ return true; }
    return false;
  }
  
  boolean contains(CellLink link){ // check if the link is contained
    return link1==link || link2==link || link3==link;
  }
  
  CellLink findLink(Cell c1, Cell c2){ // find a link between 2 cells
    for(int i=0; i < c1.links.size() ; i++){
      if(c1.links.get(i).contains(c2)){ return c1.links.get(i); }
    }
    IG.err("link not found");
    return null;
  }
  
  Cell[] oppositeCell(Cell c){ // find 2 cells opposite of a cell in a triangle
    if(cell1==c){ return  new Cell[]{ cell2, cell3 }; }
    if(cell2==c){ return  new Cell[]{ cell3, cell1 };  }
    if(cell3==c){ return  new Cell[]{ cell1, cell2 };  }
    return null;
  }
  
  CellLink oppositeLink(Cell c){ // find opposite link of a cell in a triangle
    if(!link1.contains(c)) return link1;
    if(!link2.contains(c)) return link2;
    if(!link3.contains(c)) return link3;
    return null;
  }
  
  void del(){
    cell1.faces.remove(this); // unregister this from cells and links
    cell2.faces.remove(this);
    cell3.faces.remove(this);
    link1.faces.remove(this);
    link2.faces.remove(this);
    link3.faces.remove(this);
    mesh.deleteFace(face);
  }
}