# Problem with spring force in growth algorithm

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:

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

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

``````maxForce = 100

def interact(self, agents) :
# spring force
dif = self.cell1.pos().dif(self.cell2.pos())
if force > CellLink.maxForce : # if putting too much force
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:

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

dif = self.sub(cells[l])
if force > maxForce:
force = maxForce
elif force < -maxForce:
force = -maxForce
dif.normalizeTo(force)
``````

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)
``````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
growthSpeed=0.1

IParticle.__init__(self, pos, IVec(0,0,0))
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
for a in agents :
if isinstance(a, Cell) :
if a is not self :
# push if closer than two radii
dif = a.pos().dif(self.pos())
a.push(dif)
# count neighbors and calculate their center
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

def divide(self) : # cell division
if len(self.links)==0 : # dot state
child = self.createChild(IRand.dir())
elif len(self.links)==1 : # line state
# making a triangle loop
child = self.createChild(IRand.dir())
# 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
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
dir = IVec()
dir.projectToPlane(self.nml())
child = self.createChild(dir)

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 :
return n.unit()

def createChild(self, dir) :
self.radius *= 0.5 #make both cell size half
self.pos().sub(dir) # move to the other way
self.active = False #reset activation
return child

return # no need to sort
sorted = []

def findDividingLinks(self) :  # find splittable two links and return the index numbers

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 :
removingFaces.append(f)
break
for i in range(1, len(linksToChild)-1) :
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.faces = []

def interact(self, agents) :
# spring force
dif = self.cell1.pos().dif(self.cell2.pos())
if force > CellLink.maxForce : # if putting too much force
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.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) :
# keep order of cells in right hand screw for consistent normal
self.cell1 = c2
self.cell2 = c1
self.cell3 = c3
else :
self.cell1 = c1
self.cell2 = c2
self.cell3 = c3
self.cell1.faces.append(self) # register this face to cells and links
self.cell2.faces.append(self)
self.cell3.faces.append(self)
self.face = IFace(self.cell1.vertex,self.cell2.vertex,self.cell3.vertex)

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
self.cell2 = self.cell3
self.cell3 = tmp

# 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 findLink(self, c1, c2) : # find a link between 2 cells
if l.contains(c2) :
return l
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
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)
mesh.deleteFace(self.face)
``````
my implementation (Python version / Processing 3)
``````add_library('toxiclibs')
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
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)

def draw():
background(255)

# Updating physics engine
physics.update()

for c in cells:
c.grow()
#c.interact()
c.render()

strokeWeight(1)
for i in xrange(len(physics.springs)):

# 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.id = 0
self.time = 0
self.faces = []

def grow(self):

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

# downsizing mother cell when giving birth to child cell

# location of child cell
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

# downsizing mother cell when giving birth to child cell

# location of child cell
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)

# Storing face (triangle) information in each cell
for c in (self, child, cells[self.links[0]]):

# Triangle state: forming a tetrahedron

# downsizing mother cell when giving birth to child cell

# location of child cell
child = Cell(newloc.x(), newloc.y(), newloc.z())
child.id = len(cells)

self.connect(self, child, False)

# Adding new face (triangle 1) information
for c in (self, child, cells[self.links[0]]):

# Adding new face (triangle 2) information
for c in (self, child, cells[self.links[1]]):

# Adding new face (triangle 3) information

# Tetrahedron state: edge division (core of the algorithm)
else:

# downsizing mother cell when giving birth to child cell

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

child = Cell(newloc.x(), newloc.y(), newloc.z())
child.id = len(cells)

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

def pickInterval(self):
r2 = (r1 + nInterval) % nLinks
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

if not added: # means that the child cell has not been included yet in the physics engine
cells.append(c2)

#Add spring between c1 and c2
s = VerletSpring3D(c1, c2, c2.distanceTo(c1), springForce)

#Storing linkage information for each cell

def interact(self):
maxForce = 100

dif = self.sub(cells[l])
if force > maxForce:
force = maxForce
elif force < -maxForce:
force = -maxForce
dif.normalizeTo(force)

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 )

1 Like

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.

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

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.

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.

1 Like

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

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.

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

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

nCells = 360
growthInterval = 10
divisiointervalserval = 50
repulsiveForce = .4
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)

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.id = 0
self.time = 0
self.faces = []
self.active = False

def grow(self):

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

# downsizing mother cell when giving birth to child cell

# location of child cell
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

# downsizing mother cell when giving birth to child cell

# location of child cell
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)

# Storing face (triangle) information in each cell
for c in (self, child, cells[self.links[0]]):

# Triangle state: forming a tetrahedron

# downsizing mother cell when giving birth to child cell

# location of child cell
child = Cell(newloc.x(), newloc.y(), newloc.z())
child.id = len(cells)

self.connect(self, child, False)

# Adding new face (triangle 1) information
for c in (self, child, cells[self.links[0]]):

# Adding new face (triangle 2) information
for c in (self, child, cells[self.links[1]]):

# Adding new face (triangle 3) information

# Tetrahedron state: edge division (core of the algorithm)
else:

# downsizing mother cell when giving birth to child cell

# 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
self.subSelf(dir)

child = Cell(newloc.x(), newloc.y(), newloc.z())
child.id = len(cells)

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

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

def pickInterval(self):
r2 = (r1 + nInterval) % nLinks
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

if not added: # means that the child cell has not been included yet in the physics engine
cells.append(c2)

#Add spring between c1 and c2
# s = VerletSpring3D(c1, c2, c2.distanceTo(c1), springForce)

#Storing linkage information for each cell

def interact(self):
# Spring force (Hooke's law)
dif = cells[l].sub(self)
d = dif.magnitude()
f = dif.normalize().scale(-1 * .05 * stretch)

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])
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()
``````

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

(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
``````def planarForce(self):
center, count = Vec3D(), 0