# Problem with spring force in growth algorithm

**URL:** <https://discourse.processing.org/t/problem-with-spring-force-in-growth-algorithm/10770>\
**Category:** Processing.py\
**Created:** [April 29, 2019, 11:26pm UTC](https://discourse.processing.org/t/problem-with-spring-force-in-growth-algorithm/10770 "2019-04-29T23:26:26Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [April 29, 2019, 11:26pm UTC](https://discourse.processing.org/t/problem-with-spring-force-in-growth-algorithm/10770/1 "2019-04-29T23:26:26Z")

</div>

Dear all,

I would need a hand to help me find what is wrong with the following implementation of a sketch originally relying on the iGeo library (not compatible with Processing 3).

Simply put, the algorithm consists in iteratively adding a random number of faces to a tetrahderon. The mesh gets bigger over time and the overall shape should gradually start to look like a giant 3D blob with bumps or protrusions:

 ![screenshot3](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/7/7ada252be5cd1af1acc47465fd28bfac082614e3.png)

However in my case, newly added faces tend to rapidly stand out from the mesh, forming sharp triangular spikes all over the shape.

![Capture%20d%E2%80%99%C3%A9cran%20(8)](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/2/245783df1a1a7cc2d146534622c28ceaeec3a86f.png)

see [GIF](https://media.giphy.com/media/J5vPkEpznPp9Kl7c4b/giphy.gif)

In my opinion, 2 things might cause this behavior:

- the lack of a force continuously pressuring the vertices outward (unlike my sketch, the original script uses the `push()` method from the iGeo library to pressure the vertices towards their corresponding normal)

- a difference between the physics of the `VerletSprings` class (my sketch is using the toxiclibs library) and the physics of the spring forces in the original script.

That second assumption seemed unlikely at first but I noticed that commenting out the `push()` method in the original script didn’t prevent the mesh to grow “normally” (the output looked a bit smaller in size but similar in shape, no spikes)

In order to be sure I would like to implement the snippet related to the spring force from the original script:

```auto
maxForce = 100

def interact(self, agents) : 
        # spring force
        dif = self.cell1.pos().dif(self.cell2.pos())
        force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
        if force > CellLink.maxForce : # if putting too much force
            force = CellLink.maxForce
        elif force < -CellLink.maxForce :
            force = -CellLink.maxForce
        dif.len(force)
        self.cell1.pull(dif)
        self.cell2.push(dif)

```

The problem is that I’m unable to come up with a correct transcription of that code using toxiclibs:

```auto
def interact(self):
    maxForce = 100

    for l in self.links:
        dif = self.sub(cells[l])
        force = (dif.magnitude() - (self.radius + cells[l].radius)) / (self.cell1.radius + self.cell2.radius) * 300
        if force > maxForce: 
            force = maxForce
        elif force < -maxForce:
            force = -maxForce
        dif.normalizeTo(force)
        self.addForce(dif.scale(-1))
        cells[l].addForce(dif)

```

I believe I’m calling the correct methods:

- iGeo `dif()` = toxiclibs `sub()`
- iGeo `len()` = toxiclibs `magnitude()`
- iGeo `len(force)` = toxiclibs `normalizeTo(force)` (similar to `setMag(force)` in Processing)
- I’m assuming that “pushing” or “pulling” a cell (`push()` / `pull()` in iGeo) is like adding or subtracting a force to / from that cell, so I’m using toxiclibs `addForce()`

however the result is catastrophic, almost impossible to describe (cells are violently shaking, fleeing all over the place).

**QUESTIONS** :

- Could you help me implement the `interact` function (spring force) correctly with toxiclibs ?
- Do you think something else could be the cause of that strange behavior (spikes formation) ?

* * *

> **original sketch (Python version / Processing 2 only)**
>
> ```auto
> add_library('igeo')
> 
> mesh = None # global variable of mesh geometry
> 
> def setup() : 
> size(480,360,IG.GL)
> IConfig.syncDrawAndDynamics=True #not to crash when some geometry is deleted while drawing
> mesh = IMesh().clr(0.7,0,0)
> cell = Cell(IVec(0,0,0), 10)
> 
> class Cell(IParticle) : 
> growthDuration = 1200 #duration of growth and division
> growthInterval = 10
> divisionInterval = 50
> maxRadius = 100
> growthSpeed=0.1
>     
> def __init__ (self, pos, rad) :
> IParticle. __init__ (self, pos, IVec(0,0,0))
> self.radius = rad
> self.links = [] #store links
> self.faces = [] #store faces
> self.active = False
> self.vertex = IVertex(self.pos())
> self.fric(0.2)
>   
> def interact(self, agents) : 
> neighborCenter = IVec(0,0,0)
> neighborCount=0
> neighborDist = self.radius*4 
> for a in agents : 
> if isinstance(a, Cell) : 
> if a is not self : 
> # push if closer than two radii
> if a.pos().dist(self.pos()) < a.radius+self.radius : 
> dif = a.pos().dif(self.pos())
> dif.len(((a.radius+self.radius)-dif.len())*100+50) # the closer the harder to push
> a.push(dif)
> # count neighbors and calculate their center
> if a.pos().dist(self.pos()) < a.radius+self.radius + neighborDist : 
> neighborCenter.add(a.pos())
> neighborCount+=1
> if neighborCount >= 1 : # push from center of neighbors
> neighborCenter.div(neighborCount)
> dif = self.pos().dif(neighborCenter).len(50) # constant force
> self.push(dif)
> 
> self.push(self.nml().len(20)) # pressure toward normal
>   
> def update(self) : 
> if IG.time() < Cell.growthDuration : 
> if self.time() > 0 and self.time()%Cell.divisionInterval==0 : 
> if self.active : # divide when active flag is on
> self.divide()
> if IRand.pct(50) : # random division
> self.active = True
> if self.time()%Cell.growthInterval==0 : 
> self.grow()
> if self.active : # update color by being active
> self.clr(1.0,0.5,0)
> else :
> self.clr(mesh.clr())
> self.vertex.pos().set(self.pos()) # update mesh vertex position
> self.vertex.nml(self.nml()) # update mesh vertex normal
>     
> def grow(self) : # growing cell size
> if self.radius < Cell.maxRadius : 
> self.radius += Cell.growthSpeed
>   
> def divide(self) : # cell division
> if len(self.links)==0 : # dot state
> child = self.createChild(IRand.dir())
> CellLink(self, child) # add one link
> elif len(self.links)==1 : # line state 
> # making a triangle loop
> child = self.createChild(IRand.dir())
> CellLink(child, self.links[0].cell1)
> CellLink(child, self.links[0].cell2)
> CellFace(child, self.links[0].cell1, self.links[0].cell2) # making a triangle face
> elif len(self.links)==2 : # triangle
> # making a tetrahedron enclosure
> child = self.createChild(IRand.dir())
> f = self.faces[0]
> center = IVec.center(child.pos(),f.cell1.pos(),f.cell2.pos(),f.cell3.pos())
> if f.center().dif(center).dot(f.nml())<0 :
> f.flipNml() # adjust normal to be outward
> CellLink(f.cell1, child)
> CellLink(f.cell2, child)
> CellLink(f.cell3, child)
> f1 = CellFace(f.cell1, f.cell2, child)
> f2 = CellFace(f.cell2, f.cell3, child)
> f3 = CellFace(f.cell3, f.cell1, child)
> if f1.center().dif(center).dot(f1.nml())<0 : 
> f1.flipNml()
> if f2.center().dif(center).dot(f2.nml())<0 : 
> f2.flipNml()
> if f3.center().dif(center).dot(f3.nml())<0 : 
> f3.flipNml()
> else : # link num > 3
> # pick two links and insert two faces between
> linkIdx = self.findDividingLinks()
> linksToChild = [] # links to reconnect to child
> for i in range(linkIdx[0], linkIdx[1]+1) : 
> linksToChild.append(self.links[i%len(self.links)])
> dir = IVec()
> for l in linksToChild : # average of links1 dir
> dir.add(l.oppositeDir(self).unit())
> dir.projectToPlane(self.nml())
> child = self.createChild(dir)
> self.insertChild(child, linksToChild)
>     
> def nml(self) : # calc vertex normal from face normal
> if len(self.faces)==0 :
> return IVec(0,0,1)
> n = IVec(0,0,0)
> for f in self.faces : 
> n.add(f.nml())
> return n.unit()
>     
> def createChild(self, dir) : 
> self.radius *= 0.5 #make both cell size half
> dir.len(self.radius)
> child = Cell(self.pos().cp(dir), self.radius)
> self.pos().sub(dir) # move to the other way
> self.active = False #reset activation
> return child
>     
> def sortLinks(self) : # sort links by following connected faces
> if len(self.links) <= 3 :
> return # no need to sort
> sorted = []
> currentLink = self.links[0]
> currentFace = currentLink.faces[0]
> for i in range(len(self.links)) : 
> sorted.append(currentLink)
> currentLink = currentFace.oppositeLink(currentLink.oppositeCell(self))
> currentFace = currentLink.oppositeFace(currentFace)
> self.links = sorted
>     
> def findDividingLinks(self) : # find splittable two links and return the index numbers
> linkInterval = len(self.links)//2 # interval between two links 
>         
> self.sortLinks() # sort the order of links around cell
> idx = IRand.getInt(0,len(self.links)-1)
> return [idx, idx+linkInterval] # index can be larger than self.links.size()
>     
> def insertChild(self, child, linksToChild) : # insert child cell and reconnect links
> removingFaces = []
> for f in self.faces : 
> for j in range(len(linksToChild)-1) : 
> if (f.contains(linksToChild[j]) and
> f.contains(linksToChild[j+1])) : #a face between replacing links is removed
> removingFaces.append(f)
> break
> for i in range(1, len(linksToChild)-1) : 
> linksToChild[i].delete() # replacing links are once removed
> CellLink(linksToChild[i].oppositeCell(self), child) # then recreated
> cell1 = linksToChild[0].oppositeCell(self) # one of two dividing link cell
> cell2 = linksToChild[len(linksToChild)-1].oppositeCell(self) # another dividing link cell
> CellLink(child, cell1)
> CellLink(child, cell2)
> CellLink(self, child)
> for f in removingFaces : 
> f.delete() # replace face by deleting and recreating
> cells = f.oppositeCell(self)
> CellFace(child, cells[0], cells[1])
> CellFace(self, cell1, child)
> CellFace(self, child, cell2)
> 
> class CellLink(IAgent) : # a link to connect 2 cells with spring force
> maxForce = 100
>     
> def __init__ (self, c1, c2) : 
> IAgent. __init__ (self)
> self.cell1 = c1 
> self.cell2 = c2
> self.cell1.links.append(self) # register this link to cells{
> self.cell2.links.append(self)
> self.faces = []
>     
> def interact(self, agents) : 
> # spring force
> dif = self.cell1.pos().dif(self.cell2.pos())
> force = (dif.len()-(self.cell1.radius+self.cell2.radius))/(self.cell1.radius+self.cell2.radius)*300
> if force > CellLink.maxForce : # if putting too much force
> force = CellLink.maxForce
> elif force < -CellLink.maxForce :
> force = -CellLink.maxForce
> dif.len(force)
> self.cell1.pull(dif)
> self.cell2.push(dif)
> 
> def contains(self, c) : #check if link contains the cell
> if c is self.cell1 or c is self.cell2 : 
> return True
> return False
>     
> def delete(self) : 
> self.cell1.links.remove(self) # unregister from cells
> self.cell2.links.remove(self)
> self.del() # stop agent
>   
> def oppositeFace(self, f) : # find other face on the link
> if len(self.faces)!=2 :
> return None
> if self.faces[0] is f :
> return self.faces[1]
> return self.faces[0]
>   
> def oppositeCell(self, c) : # find other cell on the link
> if self.cell1 is c : 
> return self.cell2
> if self.cell2 is c :
> return self.cell1
> print("Link does not contain the input cell")
> return None 
>   
> def oppositeDir(self, c) : #calculate a vector to the other cell
> return self.oppositeCell(c).pos().dif(c.pos())
> 
> class CellFace : # a triangle grouping 3 cells and 3 links
>     
> def __init__ (self, c1, c2, c3) : 
> self.link1 = self.findLink(c1, c2)
> # keep order of cells in right hand screw for consistent normal
> if len(self.link1.faces)==1 and self.link1.faces[0].cellOrder(c1,c2) :
> self.cell1 = c2 
> self.cell2 = c1 
> self.cell3 = c3
> else : 
> self.cell1 = c1 
> self.cell2 = c2 
> self.cell3 = c3
> self.link1 = self.findLink(self.cell1, self.cell2)
> self.link2 = self.findLink(self.cell2, self.cell3)
> self.link3 = self.findLink(self.cell3, self.cell1)
> self.cell1.faces.append(self) # register this face to cells and links
> self.cell2.faces.append(self)
> self.cell3.faces.append(self)
> self.link1.faces.append(self)
> self.link2.faces.append(self)
> self.link3.faces.append(self)
> self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)
> mesh.addFace(self.face) # add triangular face to the mesh
>     
> def center(self) : # calc center of triangle
> return IVec.center(self.cell1.pos(), self.cell2.pos(), self.cell3.pos())
>   
> def nml(self) : # normal vector
> return self.cell1.pos().nml(self.cell2.pos(), self.cell3.pos()).unit()
>     
> def flipNml(self) : # flip normal
> tmp = self.cell2
> tmpLink = self.link2
> self.cell2 = self.cell3
> self.link2 = self.link3
> self.cell3 = tmp
> self.link3 = tmpLink
> 
> # True if ordr of c1-c2 follows the order of cell1-cell2-cell3, otherwise False
> def cellOrder(self, c1, c2) : 
> if ( self.cell1 is c1 and self.cell2 is c2 or 
> self.cell2 is c1 and self.cell3 is c2 or 
> self.cell3 is c1 and self.cell1 is c2 ) : 
> return True
> return False
>     
> def contains(self, link) : # check if the link is contained
> return self.link1 is link or self.link2 is link or self.link3 is link
>     
> def findLink(self, c1, c2) : # find a link between 2 cells
> for l in c1.links : 
> if l.contains(c2) : 
> return l
> print("link not found")
> return None
>   
> def oppositeCell(self, c) : # find 2 cells opposite of a cell in a triangle
> if self.cell1 is c :
> return [self.cell2, self.cell3]
> if self.cell2 is c :
> return [self.cell3, self.cell1]
> if self.cell3 is c :
> return [self.cell1, self.cell2]
> return None
>   
> def oppositeLink(self, c) : # find opposite link of a cell in a triangle
> if not self.link1.contains(c) :
> return self.link1
> if not self.link2.contains(c) :
> return self.link2
> if not self.link3.contains(c) :
> return self.link3
> return None
>     
> def delete(self) : 
> self.cell1.faces.remove(self) # unregister self from cells and links
> self.cell2.faces.remove(self)
> self.cell3.faces.remove(self)
> self.link1.faces.remove(self)
> self.link2.faces.remove(self)
> self.link3.faces.remove(self)
> mesh.deleteFace(self.face)
> 
> ```

> **my implementation (Python version / Processing 3)**
>
> ```auto
> add_library('toxiclibs')
> add_library('peasycam')
> from toxi.physics3d.behaviors import GravityBehavior3D, AttractionBehavior3D
> from toxi.physics3d import VerletPhysics3D, VerletParticle3D, VerletSpring3D
> import random as rdm
> 
> nCells = 360
> growthInterval = 10
> divisiointervalserval = 50
> repulsiveForce = .4
> nRadii = 5
> springForce = .001
> 
> colors = [(random(255), random(255), random(255)) for i in xrange(nCells)]
> 
> def setup():
> perspective(60 * DEG_TO_RAD, width/float(height), 2, 6000)
> size(700, 400, P3D)
> fill(255, 40, 70)
> smooth(8)
>     
> global physics, cells
>     
> # Array containing the "Cell" objects
> cells = [Cell(0, 0, 0)]
>     
> cam = PeasyCam(this, 80)
> cam.setWheelScale(.3)
>     
> physics = VerletPhysics3D()
> physics.setDrag(.1)
> physics.addParticle(cells[0])
> physics.addBehavior(AttractionBehavior3D(cells[0], cells[0].radius*nRadii, -repulsiveForce))
>     
>     
> def draw():
> background(255)
>     
> # Updating physics engine
> physics.update()
> 
> for c in cells:
> c.grow()
> #c.interact()
> c.render()
>         
>         
> # Adjusting spring length (sum of the radii of the 2 cells forming a link)
> strokeWeight(1)
> for i in xrange(len(physics.springs)):
> physics.springs.get(i).setRestLength(physics.springs.get(i).a.radius + physics.springs.get(i).b.radius)
>         
> # Drawing springs
> line(physics.springs.get(i).a.x(), physics.springs.get(i).a.y(), physics.springs.get(i).a.z(), physics.springs.get(i).b.x(), physics.springs.get(i).b.y(), physics.springs.get(i).b.z())
>     
>     
> class Cell(VerletParticle3D):    
> def __init__ (self, x, y, z):
> super(Cell, self). __init__ (x, y, z)
> self.radius = 5
> self.maxRadius = 20
> self.id = 0
> self.time = 0
> self.links = []
> self.faces = []
>         
> def grow(self):
>         
> # Incrementing cell's lifetime
> self.time += 1
>         
> # As long as we haven't reached the max number of cells yet --> keep growing
> if len(cells) < nCells:
> if random(1) > .5:
> if self.time > 0 and self.time%divisiointervalserval == 0:
>                     
> # Point state: forming a line
> if not self.links:
>                         
> # downsizing mother cell when giving birth to child cell
> self.radius *= .5
>                         
> # location of child cell
> newloc = self.add(Vec3D.randomVector().scale(self.radius))
> child = Cell(newloc.x(), newloc.y(), newloc.z())
> child.id = len(cells)
>                         
> # creating a link between mother cell and child cell
> self.connect(self, child, False)
>                      
>                          
> # Line state: forming a triangle
> elif len(self.links) == 1:
>                                                 
> # downsizing mother cell when giving birth to child cell
> self.radius *= .5
>                         
> # location of child cell
> newloc = self.add(Vec3D.randomVector().normalizeTo(self.radius))
> child = Cell(newloc.x(), newloc.y(), newloc.z())
> child.id = len(cells)
>                         
> # creating a link between mother cell and child cell
> # creating a link between child cell and cell on the other side of the link
> self.connect(self, child, False)
> self.connect(child, cells[self.links[0]], True)
>                         
> # Storing face (triangle) information in each cell
> for c in (self, child, cells[self.links[0]]):
> c.faces.append(frozenset([self.id, child.id, cells[self.links[0]].id]))
>                      
>                    
> # Triangle state: forming a tetrahedron  
> elif len(self.links) == 2:
>                         
> # downsizing mother cell when giving birth to child cell
> self.radius *= .5
>                         
> # location of child cell
> newloc = self.computeNormal([self.id, cells[self.links[0]].id, cells[self.links[1]].id])
> child = Cell(newloc.x(), newloc.y(), newloc.z())
> child.id = len(cells)
>                         
> self.connect(self, child, False)
> self.connect(cells[self.links[0]], child, True)
> self.connect(cells[self.links[1]], child, True)
>                         
> # Adding new face (triangle 1) information
> for c in (self, child, cells[self.links[0]]):
> c.faces.append(frozenset([self.id, child.id, cells[self.links[0]].id]))
>                             
> # Adding new face (triangle 2) information
> for c in (self, child, cells[self.links[1]]):
> c.faces.append(frozenset([self.id, child.id, cells[self.links[1]].id]))
>                         
> # Adding new face (triangle 3) information
> for c in (cells[self.links[0]], child, cells[self.links[1]]):
> c.faces.append(frozenset([cells[self.links[0]].id, child.id, cells[self.links[1]].id]))
>                         
>                         
> # Tetrahedron state: edge division (core of the algorithm)
> else:
>                         
> # downsizing mother cell when giving birth to child cell
> self.radius *= .5
>                         
> # Picking 2 links at random and computing the interval between (number of faces)
> intervals, vertices = self.pickInterval()
> print intervals
>                         
> # Computing child cell location
> n = self.computeNormal(list(vertices)) # centered normal pointing outward
> dir = n.sub(self)
> newloc = self.add(dir.normalizeTo(self.radius))
>                     
> child = Cell(newloc.x(), newloc.y(), newloc.z())
> child.id = len(cells)
>                         
> # Creating new links 
> for i, v in enumerate(vertices):
> if i == 0: self.connect(self, child, False)
> else: self.connect(cells[v], child, True)
>                                 
> # Adding new face (triangle) information
> for c in (cells[vertices[i]], child, cells[vertices[(i+1)%len(vertices)]]):
> c.faces.append(frozenset([cells[vertices[i]].id, child.id, cells[vertices[(i+1)%len(vertices)]].id]))
>                           
>                         
> # Remove former links (info + corresponding spring) between self and its selected neighboring cell
> if intervals > 1:      
> for v in vertices[2: -1]:
> self.links.remove(cells[v].id)
> cells[v].links.remove(self.id)
> physics.removeSpring(physics.getSpring(self, cells[v]))
>                                 
>                                 
> # Removing former face (triangle) information from each cell
> for i, v in enumerate(vertices):
> if i > 0 and i < len(vertices) - 1:
> for c in (self, cells[vertices[i]], cells[vertices[i+1]]):
> c.faces.remove(frozenset([self.id, cells[vertices[i]].id, cells[vertices[i+1]].id]))      
>                           
> 
> # Increasing cell size intermittently
> if self.time%growthInterval == 0:
> if self.radius < self.maxRadius:
> self.radius += .1
>                     
>     
> def pickInterval(self):
> nLinks = len(self.links)
> nInterval = nLinks // 2
> r1 = int(random(nLinks))
> r2 = (r1 + nInterval) % nLinks
> l1 = self.links[r1]
> l2 = self.links[r2]
> vertices = [self.id, l1]
> oppcell = l1
> visited = []
>             
> while oppcell != l2:
> faces = [f for f in cells[oppcell].faces if self.id in f]
> sface = faces[0] if faces[0] not in visited else faces[1]
> visited.append(sface)
> oppcell = iter(sface - set((self.id, oppcell))).next()
> vertices.append(oppcell)
>     
> interval = len(vertices) - 2
> 
> return interval, vertices                    
>     
>     
>     
> def connect(self, c1, c2, added):
>         
> if not added: # means that the child cell has not been included yet in the physics engine
> cells.append(c2)
> physics.addParticle(c2)
> physics.addBehavior(AttractionBehavior3D(c2, c2.radius*nRadii, -repulsiveForce))
>         
> #Add spring between c1 and c2
> s = VerletSpring3D(c1, c2, c2.distanceTo(c1), springForce) 
> physics.addSpring(s)  
>         
> #Storing linkage information for each cell
> c1.links.append(c2.id)
> c2.links.append(c1.id)  
>         
> def interact(self):
> maxForce = 100
> 
> for l in self.links:
> dif = self.sub(cells[l])
> force = (dif.magnitude() - (self.radius + cells[l].radius)) / (self.radius + cells[l].radius) * 300
> if force > maxForce: 
> force = maxForce
> elif force < -maxForce:
> force = -maxForce
> dif.normalizeTo(force)
> self.addForce(dif.scale(-1))
> cells[l].addForce(dif)
>             
>         
> def computeNormal(self, indices):
>         
> # Calculate normal with Newell's method
> nml = PVector()
>         
> for i in xrange(len(indices)):
>                     
> p0 = cells[indices[i]] #current vertex
> p1 = cells[indices[(i+1)%len(indices)]] #next vertex
>             
> nml.x += (p0.y() - p1.y()) * (p0.z() + p1.z())
> nml.y += (p0.z() - p1.z()) * (p0.x() + p1.x())
> nml.z += (p0.x() - p1.x()) * (p0.y() + p1.y())
>             
> # Normalizing then scaling-up to see a line
> nml = nml.normalize()#.mult(1.1)
>             
> # Center of the polygon
> center = PVector()
> for i in range(len(indices)):
> center.x += cells[indices[i]].x()
> center.y += cells[indices[i]].y()
> center.z += cells[indices[i]].z()
>                 
> center.div(len(indices))
>         
> # Centroid of tetrahedron
> centroid = PVector()
>             
> # Check if normal is heading inwards or outwards
> d = PVector.dot(center - centroid, nml)
>         
> if d < 0: nml.mult(-1)
>              
> return Vec3D(center.x + nml.x, center.y + nml.y, center.z + nml.z)
> 
>         
> def render(self):
> strokeWeight(4)
> point(self.x(), self.y(), self.z())
>     
> beginShape(TRIANGLE)
> strokeWeight(1)
> stroke(0)
> for f in self.faces:
> f = list(f)
> fill(colors[f[0]][0], colors[f[1]][1], colors[f[2]][2])
> vertex(cells[f[0]].x(), cells[f[0]].y(), cells[f[0]].z())
> vertex(cells[f[1]].x(), cells[f[1]].y(), cells[f[1]].z())
> vertex(cells[f[2]].x(), cells[f[2]].y(), cells[f[2]].z())
> endShape()
> 
> ```

(if you want to try the `interact` function (line 44), line 219 and 51 need to be commented out to avoid calling the `VerletSpring` class at the same time )

---

<div class="post-metadata">

**Author:** ![jeremydouglass](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/jeremydouglass/32/20_2.png) [@jeremydouglass](https://discourse.processing.org/u/jeremydouglass)\
**Post date:** [April 30, 2019, 5:15am UTC](https://discourse.processing.org/t/problem-with-spring-force-in-growth-algorithm/10770/2 "2019-04-30T05:15:18Z")

</div>

From visual inspection, a guess: your upper shape “grows” by splitting two joined faces into four faces by adding a center point to their shared edge, then extruding that point for all four faces.

Yet it looks like your lower image isn’t geometrically the same. It appears to have geometry like skin tags – a group of faces held on to the main body by a single shared edge. If you were using a divide algorithm like above, that should be impossible.

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [April 30, 2019, 5:26pm UTC](https://discourse.processing.org/t/problem-with-spring-force-in-growth-algorithm/10770/3 "2019-04-30T17:26:23Z")

</div>

Hi @jeremydouglass, thank you for the reply.

Your first guess is very close to the actual algorithm, the mesh “grows” by:

- selecting 2 edges on a dividing cell
- inserting a new child cell
- creating 3 new edges and 2 new faces

![](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/6/6518570dc643e4346e3e18df5f7ffdf4af37ec09.png)  
 ![](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/6/603bdd6159f59e35410b256520ec7f8f19a7d521.png)  
 ![](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/6/68528851fb1b2bb7ed42e9375d3009248b46bfdd.png)

> [@jeremydouglass](#):
>
> Yet it looks like your lower image isn’t geometrically the same […] If you were using a divide algorithm like above, that should be impossible.

I have checked, double-checked and triple-checked and I am confident that the implementation of the dividing process _is_ correct. If there’s a mistake (and there must be one!) it should be somewhere else in the code.

In the original sketch it looks like the newly added cells **stay on the surface of the shape** whereas in my case they immediatly stand out from the mesh.  
I’m pretty sure that doesn’t have to do with the repulsive force (cells have a repulsive force to prevent overlaps), neither it has to do with the birth location.

My take is that the physics of the spring connections is critical in the shaping of the mesh but I’m still unable to implement the formula used in the original sketch.

---

<div class="post-metadata">

**Author:** ![jb4x](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/jb4x/32/789_2.png) [@jb4x](https://discourse.processing.org/u/jb4x)\
**Post date:** [May 1, 2019, 4:21pm UTC](https://discourse.processing.org/t/problem-with-spring-force-in-growth-algorithm/10770/4 "2019-05-01T16:21:25Z")

</div>

Hi @solub,

I tried to mess around with it but I don’t know ToxicLib and never used Python before so you can imagine I did not have lots of success ^^

But, I may have some leads for you (maybe):

- When you create a new cell and compute the new location, I added +5 in the scale factor so the new particule is not overlapping the old one.
- I displayed each particles as sphere instead of point and deleted the triangles to have a better view of what was going on and I noticed that when the particles are at a resting place (I limited the number of cells to be 2) they are not touching each others. So there is probably something wrong with the way the springs are working.
- Something that I don’t understand is that you add some string but you also add some behavior. Shouldn’t be the strings enough ?
- I locked the initial particle so it will never move because when playing with the strings and behavior the particles started to fly all over the place…

Also, one thing that annoyed me is that it seems that the simulation and the display are linked together. I wanted to create a mousePressed function in which I could add a new particle anytime I want and just keep the simulation running in the draw loop but I couldn’t do it.

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [May 1, 2019, 5:46pm UTC](https://discourse.processing.org/t/problem-with-spring-force-in-growth-algorithm/10770/5 "2019-05-01T17:46:07Z")

</div>

– @jb4x thanks for the reply, I really appreciate you took the time to help.

> [@jb4x](#):
>
> Something that I don’t understand is that you add some string but you also add some behavior. Shouldn’t be the strings enough ?

Indeed each particle is given a repulsive force to prevent collisions / overlapping (that’s why they’re not touching). On the other hand we don’t want the particles to get too far from each others so we connect them with springs.

> [@jb4x](#):
>
> I wanted to create a mousePressed function in which I could add a new particle anytime I want…

Here’s a version with spheres and a `mousePressed` function (just replaced line 68 but it’s not ideal, a `mouseReleased` version should be much more convenient).

```auto
add_library('toxiclibs')
add_library('peasycam')
from toxi.physics3d.behaviors import GravityBehavior3D, AttractionBehavior3D
from toxi.physics3d import VerletPhysics3D, VerletParticle3D, VerletSpring3D

nCells = 360
growthInterval = 10
divisiointervalserval = 50
repulsiveForce = .4
nRadii = 5
springForce = .001

colors = [(random(255), random(255), random(255)) for i in xrange(nCells)]

def setup():
    perspective(60 * DEG_TO_RAD, width/float(height), 2, 6000)
    sphereDetail(9)
    size(700, 400, P3D)
    #frameRate(3)
    fill(255, 40, 70)
    smooth(8)
    
    global physics, cells
    
    # Array containing the "Cell" objects
    cells = [Cell(0, 0, 0)]
    
    cam = PeasyCam(this, 80)
    cam.setWheelScale(.1)
    
    physics = VerletPhysics3D()
    physics.setDrag(.4)
    physics.addParticle(cells[0])
    physics.addBehavior(AttractionBehavior3D(cells[0], cells[0].radius*nRadii, -repulsiveForce))
    
    
def draw():
    background(255)
    
    
    # Updating physics engine
    physics.update()
    
    for c in cells:
        c.grow()
        c.interact()
        c.render()
        

class Cell(VerletParticle3D):    
    def __init__ (self, x, y, z):
        super(Cell, self). __init__ (x, y, z)
        self.radius = 5
        self.maxRadius = 25
        self.id = 0
        self.time = 0
        self.links = []
        self.faces = []
        self.active = False
        
    def grow(self):
        
        # Incrementing cell's lifetime
        self.time += 1
        
        # As long as we haven't reached the max number of cells yet --> keep growing
        if len(cells) < nCells:
            if keyPressed:
                if self.active:
                    
                    # Point state: forming a line
                    if not self.links:
                        
                        # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        
                        # location of child cell
                        newloc = self.add(Vec3D.randomVector().scale(self.radius))
                        child = Cell(newloc.x(), newloc.y(), newloc.z())
                        child.id = len(cells)
                        
                        # creating a link between mother cell and child cell
                        self.connect(self, child, False)
                     
                         
                    # Line state: forming a triangle
                    elif len(self.links) == 1:
                                                
                        # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        
                        # location of child cell
                        newloc = self.add(Vec3D.randomVector().normalizeTo(self.radius))
                        child = Cell(newloc.x(), newloc.y(), newloc.z())
                        child.id = len(cells)
                        
                        # creating a link between mother cell and child cell
                        # creating a link between child cell and cell on the other side of the link
                        self.connect(self, child, False)
                        self.connect(child, cells[self.links[0]], True)
                        
                        # Storing face (triangle) information in each cell
                        for c in (self, child, cells[self.links[0]]):
                            c.faces.append(frozenset([self.id, child.id, cells[self.links[0]].id]))
                     
                   
                    # Triangle state: forming a tetrahedron  
                    elif len(self.links) == 2:
                        
                         # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        
                        # location of child cell
                        newloc = self.computeNormal([self.id, cells[self.links[0]].id, cells[self.links[1]].id])
                        child = Cell(newloc.x(), newloc.y(), newloc.z())
                        child.id = len(cells)
                        
                        self.connect(self, child, False)
                        self.connect(cells[self.links[0]], child, True)
                        self.connect(cells[self.links[1]], child, True)
                        
                        # Adding new face (triangle 1) information
                        for c in (self, child, cells[self.links[0]]):
                            c.faces.append(frozenset([self.id, child.id, cells[self.links[0]].id]))
                            
                        # Adding new face (triangle 2) information
                        for c in (self, child, cells[self.links[1]]):
                            c.faces.append(frozenset([self.id, child.id, cells[self.links[1]].id]))
                        
                        # Adding new face (triangle 3) information
                        for c in (cells[self.links[0]], child, cells[self.links[1]]):
                            c.faces.append(frozenset([cells[self.links[0]].id, child.id, cells[self.links[1]].id]))
                        
                        
                    # Tetrahedron state: edge division (core of the algorithm)
                    else:
                        
                        # downsizing mother cell when giving birth to child cell
                        self.radius *= .5
                        
                        # Picking 2 links at random and computing the interval between (number of faces)
                        intervals, vertices = self.pickInterval()
                                                
                        # Computing child cell location
                        n = self.computeNormal(list(vertices)) # centered normal pointing outward
                        dir = n.normalizeTo(self.radius)
                        self.subSelf(dir)
                        newloc = self.add(dir)
                    
                        child = Cell(newloc.x(), newloc.y(), newloc.z())
                        child.id = len(cells)
                        
                        
                        # Creating new links 
                        for i, v in enumerate(vertices):
                            if i == 0: self.connect(self, child, False)
                            else: self.connect(cells[v], child, True)
                                
                            # Adding new face (triangle) information
                            for c in (cells[vertices[i]], child, cells[vertices[(i+1)%len(vertices)]]):
                                c.faces.append(frozenset([cells[vertices[i]].id, child.id, cells[vertices[(i+1)%len(vertices)]].id]))
                          
                        
                        # Remove former links (info + corresponding spring) between self and its selected neighboring cell
                        if intervals > 1:      
                            for v in vertices[2: -1]:
                                self.links.remove(cells[v].id)
                                cells[v].links.remove(self.id)
                                
                                
                        # Removing former face (triangle) information from each cell
                        for i, v in enumerate(vertices):
                            if i > 0 and i < len(vertices) - 1:
                                for c in (self, cells[vertices[i]], cells[vertices[i+1]]):
                                    c.faces.remove(frozenset([self.id, cells[vertices[i]].id, cells[vertices[i+1]].id]))      
                          
                        self.active = False
                        
                        
                if random(1) > .5:
                    self.active = True
                    
                    
                    
            # Increasing cell size intermittently
            if self.time%growthInterval == 0:
                if self.radius < self.maxRadius:
                    self.radius += .1
                                        
    
    def pickInterval(self):
        nLinks = len(self.links)
        nInterval = nLinks // 2
        r1 = int(random(nLinks))
        r2 = (r1 + nInterval) % nLinks
        l1 = self.links[r1]
        l2 = self.links[r2]
        vertices = [self.id, l1]
        oppcell = l1
        visited = []
            
        while oppcell != l2:
            faces = [f for f in cells[oppcell].faces if self.id in f]
            sface = faces[0] if faces[0] not in visited else faces[1]
            visited.append(sface)
            oppcell = iter(sface - set((self.id, oppcell))).next()
            vertices.append(oppcell)
    
        interval = len(vertices) - 2

        return interval, vertices                    
    
    
    
    def connect(self, c1, c2, added):
        
        if not added: # means that the child cell has not been included yet in the physics engine
            cells.append(c2)
            physics.addParticle(c2)
            physics.addBehavior(AttractionBehavior3D(c2, c2.radius*nRadii, -repulsiveForce))
        
        #Add spring between c1 and c2
        # s = VerletSpring3D(c1, c2, c2.distanceTo(c1), springForce) 
        # physics.addSpring(s)  
        
        #Storing linkage information for each cell
        c1.links.append(c2.id)
        c2.links.append(c1.id)  
        
    def interact(self):
        # Spring force (Hooke's law)
        for l in self.links:
            dif = cells[l].sub(self) 
            d = dif.magnitude()
            stretch = (d - (self.radius + cells[l].radius)) #/ (self.radius + cells[l].radius) * 300
            f = dif.normalize().scale(-1 * .05 * stretch)
            cells[l].addForce(f)

        
    def computeNormal(self, indices):
        
        # Calculate normal with Newell's method
        nml = PVector()
        
        for i in xrange(len(indices)):
                    
            p0 = cells[indices[i]] #current vertex
            p1 = cells[indices[(i+1)%len(indices)]] #next vertex
            
            nml.x += (p0.y() - p1.y()) * (p0.z() + p1.z())
            nml.y += (p0.z() - p1.z()) * (p0.x() + p1.x())
            nml.z += (p0.x() - p1.x()) * (p0.y() + p1.y())
            
        # Normalizing then scaling-up to see a line
        nml = nml.normalize()#.mult(1.1)
            
        # Center of the polygon
        center = PVector()
        for i in range(len(indices)):
            center.x += cells[indices[i]].x()
            center.y += cells[indices[i]].y()
            center.z += cells[indices[i]].z()
                
        center.div(len(indices))
        
        # Centroid of tetrahedron
        centroid = PVector()
            
        # Check if normal is heading inwards or outwards
        d = PVector.dot(center - centroid, nml)
        
        if d < 0: nml.mult(-1)
             
        #return Vec3D(center.x, center.y, center.z)
        return Vec3D(center.x + nml.x, center.y + nml.y, center.z + nml.z)

        
    def render(self):
        strokeWeight(4)
        point(self.x(), self.y(), self.z())
        
        pushMatrix()
        translate(self.x(), self.y(), self.z())
        strokeWeight(1)
        noStroke()
        fill(colors[self.id][0], colors[self.id][1], colors[self.id][2])
        sphere(self.radius)
        popMatrix()
    
        beginShape(TRIANGLE)
        strokeWeight(1)
        stroke(0)
        noFill()
        for f in self.faces:
            f = list(f)
            vertex(cells[f[0]].x(), cells[f[0]].y(), cells[f[0]].z())
            vertex(cells[f[1]].x(), cells[f[1]].y(), cells[f[1]].z())
            vertex(cells[f[2]].x(), cells[f[2]].y(), cells[f[2]].z())
        endShape()

```

![Capture%20d%E2%80%99%C3%A9cran%20(20)](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/4/4b6d3e8b1d676e29fcaa90153507f1c0e032c9f0.png)

Also I’ve replaced toxiclibs `VerletSprings` by an `interact` function that computes the same Hooke’s law. The result is the same obviously but what I’d like to do now is implement the original sketch springs calculations (a bit different from Hooke’s law) (see the original question).

---

<div class="post-metadata">

**Author:** ![solub](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.processing.org/solub/32/333_2.png) [@solub](https://discourse.processing.org/u/solub)\
**Post date:** [May 11, 2019, 12:39pm UTC](https://discourse.processing.org/t/problem-with-spring-force-in-growth-algorithm/10770/6 "2019-05-11T12:39:32Z")

</div>

_(Just realized I forgot to update the thread with the solution I had found)_

The problem was towfold:

- I had to keep the first and last edges instead of removing/rebuilding them (connecting order is important)
- Planar force was missing

```auto
def planarForce(self):
    center, count = Vec3D(), 0
    for l in self.links:
        center.addSelf(cells[l])
        count += 1
                
    center.scaleSelf(1.0/count)
    dif = center.sub(self)
    f = dif.sub(self.getVelocity().scale(.3))
    self.addForce(f)

```

 ![screenshot%20(2)](https://canada1.discourse-cdn.com/flex036/uploads/processingfoundation1/original/2X/1/116029572c8b74e94fe3110ff4d59b9bea362bed.jpeg)  
[gif](https://media.giphy.com/media/QyJ1sJsQSWZFD8177f/giphy.gif)