Block Research Group BRG code eQUILIBRIUM COMPAS

A0. Force Density Method

Form Finding Principle

Form Finding with the Force Density Method

For reference, the force density method is based on

Schek H., The Force Density Method for Form Finding and Computation of General Networks, Computer Methods in Applied Mechanics and Engineering 3: 115-134, 1974.

and nicely explained in

Linkwitz K., Force density method - Design of a timber shell, Shell Structures for Architecture: Form Finding and Optimization, Adriaenssens, S., Block, P., Veenendaal, D. and Williams, C. (editors), Routledge, London, 2014.

However, we will go through the principle with its basic steps here:

The Steps

a. network

First, we will create a network data structure with a random geometry that is far from its equilibrium state:

  • The network will have default attributes for its nodes and edges, such as the information if a node is an anchor, the nodal residual forces in all x-, y-, and z-direction, and the force in the edges. These will be set as default to

    • is_anchor=False, meaning not an anchor,

    • rx,yz=0, meaning no residual forces yet, and

    • f=1, meaning uniform edge force.

  • The network will consist of five nodes (a to e) of which the leaves are anchored, and of edges that all connect to the middle node.

  • The network is visualised with its anchors in red. 

from compas.datastructures import Network

import compas_rhino
from compas_rhino.artists import NetworkArtist

# ==============================================================================
# create a network
# ==============================================================================

network = Network()

network.update_default_node_attributes(is_anchor=False)
network.update_default_node_attributes(rx=0, ry=0, rz=0)
network.update_default_edge_attributes(f=1)

a = network.add_node(x=0, y=0, z=0, is_anchor=True)
b = network.add_node(x=10, y=0, z=10, is_anchor=True)
c = network.add_node(x=10, y=10, z=0, is_anchor=True)
d = network.add_node(x=0, y=10, z=10, is_anchor=True)

e = network.add_node(x=5, y=5, z=0)

network.add_edge(a, e)
network.add_edge(b, e)
network.add_edge(c, e)
network.add_edge(d, e)

# ==============================================================================
# visualize the geometry
# ==============================================================================

# first clear the Rhino model
compas_rhino.clear()

artist = NetworkArtist(network)
artist.layer = "CSD2::A0::FormFinding"

# color the anchors red
node_color = {node: (255, 0, 0) for node in network.nodes_where({'is_anchor': True})}

artist.draw_nodes(color=node_color)
artist.draw_edges()

b. forces

The edges pull on the nodes with a force of 1 unit, as defined in the default edge attributes. These forces are visualised by red arrows. 

# ==============================================================================
# visualize the forces 
# ==============================================================================

lines = []
for node in network.nodes():
    a = network.node_attributes(node, 'xyz')

    for nbr in network.neighbors(node):
        b = network.node_attributes(nbr, 'xyz')

        edge = node, nbr
        if not network.has_edge(*edge):
            edge = nbr, node

        force = network.edge_attribute(edge, 'f')
        direction = normalize_vector(subtract_vectors(b, a))
        vector = scale_vector(direction, force)

        lines.append({
            'start': a,
            'end': add_vectors(a, vector),
            'arrow': 'end',
            'color': (255, 0, 0)
        })

compas_rhino.draw_lines(lines, layer=artist.layer)

c. residuals

The sum of all forces in the edges (and external loads) acting on a node are the residual forces; it is basically the leftover force in the free nodes (shown in cyan). The residual forces on the anchors are the reaction forces (shown in green).

If all these edge forces acting on one free node were in equilibrium, the sum of these forces would be equal to zero in all three x-, y- and z-direction:

F1,x+F2,x+F3,x+F4,x=0F1,y+F2,y+F3,y+F4,y=0F1,z+F2,z+F3,z+F4,z=0F_{1, x} + F_{2, x} + F_{3, x} + F_{4, x} = 0 \\F_{1, y} + F_{2, y} + F_{3, y} + F_{4, y} = 0 \\F_{1, z} + F_{2, z} + F_{3, z} + F_{4, z} = 0

However, for a random geometry, this is very unlikely, and the edge forces would sum up to a residual force R (in red):

F1,x+F2,x+F3,x+F4,x=RxF1,y+F2,y+F3,y+F4,y=RyF1,z+F2,z+F3,z+F4,z=RzF_{1, x} + F_{2, x} + F_{3, x} + F_{4, x} = R_x \\F_{1, y} + F_{2, y} + F_{3, y} + F_{4, y} = R_y \\F_{1, z} + F_{2, z} + F_{3, z} + F_{4, z} = R_z

This can also be expressed in a more compact manner with the sum symbol:

∑j=14Fj,x=Rx∑j=14Fj,y=Ry∑j=14Fj,z=Rz\sum_{j=1}^4 F_{j, x} = R_x \\\sum_{j=1}^4 F_{j, y} = R_y \\\sum_{j=1}^4 F_{j, z} = R_z

The force components in the x-, y- and z-direction can be expressed as follows: (e.g. for the x-component the projection to the x-axis is the difference in x-coordinates divided by the edge length.)

∑j=14Fjxj−x0Lj=Rx∑j=14Fjyj−y0Lj=Ry∑j=14Fjzj−z0Lj=Rz\sum_{j=1}^4 F_{j}\frac{x_j - x_0}{L_j} = R_x \\\sum_{j=1}^4 F_{j}\frac{y_j - y_0}{L_j} = R_y \\\sum_{j=1}^4 F_{j}\frac{z_j - z_0}{L_j} = R_z

So let's compute the residual forces: We define the residual as a force vector that is zero. We then add to it the residual forces from each of the neighbouring edges. Computing the forces in the edges is the same as in b.forces, but then we must separate the forces into its x-, y-, and z-components. We store the sum of residuals to the nodes by updating its attributes.

# ==============================================================================
# compute the residual forces in the current geometry
# ==============================================================================

for node in network.nodes():
    A = network.node_attributes(node, 'xyz')
    R = [0, 0, 0]
    for nbr in network.neighbors(node):
        B = network.node_attributes(nbr, 'xyz')
        edge = node, nbr
        if not network.has_edge(*edge):
            edge = nbr, node
        F = network.edge_attribute(edge, 'f')
        L = network.edge_length(*edge)
        R[0] += F * (B[0] - A[0]) / L
        R[1] += F * (B[1] - A[1]) / L
        R[2] += F * (B[2] - A[2]) / L
    network.node_attributes(node, ['rx', 'ry', 'rz'], R)

And let's also visualize the residual forces. The residuals at the anchor points are the reaction forces and are thus displayed in green and with flipped orientation: 

# ==============================================================================
# visualize the reaction forces >> NEW (no more forces)
# ==============================================================================

lines = []
for node in network.nodes_where({'is_anchor': True}):
    start = network.node_attributes(node, 'xyz')
    residual = network.node_attributes(node, ['rx', 'ry', 'rz'])
    end = subtract_vectors(start, residual)
    color = (0, 255, 0)
    lines.append({'start': start, 'end': end, 'arrow': 'end', 'color': color})

compas_rhino.draw_lines(lines, layer=layer)

# ==============================================================================
# visualize the residual forces >> NEW (no more forces)
# ==============================================================================

lines = []
for node in network.nodes_where({'is_anchor': False}):
    start = network.node_attributes(node, 'xyz')
    residual = network.node_attributes(node, ['rx', 'ry', 'rz'])
    end = add_vectors(start, residual)
    color = (0, 255, 255)
    lines.append({'start': start, 'end': end, 'arrow': 'end', 'color': color})

compas_rhino.draw_lines(lines, layer=layer)

d. update geometry

So if the residuals are not zero (or below a threshold), the node must move towards an equilibrium configuration. So we start moving the node in direction of the residual force, and by doing that the residual force becomes smaller. In order not to shoot too far, we will add a damping factor of 0.5.

Remember to recompute the new residual force after updating the geometry as it will have changed!

# ==============================================================================
# compute the residual forces in the current geometry
# ==============================================================================
# ...before
# ==============================================================================
# move free nodes in direction of residual > NEW!!!
# ==============================================================================

for node in network.nodes():
    if network.node_attribute(node, 'is_anchor'):
        continue

    rx, ry, rz = network.node_attributes(node, ['rx', 'ry', 'rz'])
    x0, y0, z0 = network.node_attributes(node, 'xyz')
    x1 = x0 + 0.5 * rx
    y1 = y0 + 0.5 * ry
    z1 = z0 + 0.5 * rz
    network.node_attributes(node, 'xyz', [x1, y1, z1])

# ...afterwards
# ==============================================================================
# recompute all residuals in the new geometry > NEW TO BE REPEATED!!!
# ==============================================================================

for node in network.nodes():
    A = network.node_attributes(node, 'xyz')
    R = [0, 0, 0]
    for nbr in network.neighbors(node):
        B = network.node_attributes(nbr, 'xyz')
        edge = node, nbr
        if not network.has_edge(*edge):
            edge = nbr, node
        F = network.edge_attribute(edge, 'f')
        L = network.edge_length(*edge)
        R[0] += F * (B[0] - A[0]) / L
        R[1] += F * (B[1] - A[1]) / L
        R[2] += F * (B[2] - A[2]) / L
    network.node_attributes(node, ['rx', 'ry', 'rz'], R)

e. equilibrium procedure

We must repeat the previous steps until equilibrium is found, meaning the residuals are zero (or below a threshold).

In order to keep our script clean and not to repeat code, we first pack the code blocks into functions. 

# ==============================================================================
# helpers
# ==============================================================================

def update_residuals(network):
    for node in network.nodes():
        A = network.node_attributes(node, 'xyz')
        R = [0, 0, 0]
        for nbr in network.neighbors(node):
            B = network.node_attributes(nbr, 'xyz')
            edge = node, nbr
            if not network.has_edge(*edge):
                edge = nbr, node
            F = network.edge_attribute(edge, 'f')
            L = network.edge_length(*edge)
            R[0] += F * (B[0] - A[0]) / L
            R[1] += F * (B[1] - A[1]) / L
            R[2] += F * (B[2] - A[2]) / L
        network.node_attributes(node, ['rx', 'ry', 'rz'], R)


def update_geometry(network):
    for node in network.nodes_where({'is_anchor': False}):
        rx, ry, rz = network.node_attributes(node, ['rx', 'ry', 'rz'])
        x0, y0, z0 = network.node_attributes(node, 'xyz')
        x1 = x0 + 0.5 * rx
        y1 = y0 + 0.5 * ry
        z1 = z0 + 0.5 * rz
        network.node_attributes(node, 'xyz', [x1, y1, z1])


def draw_reactions(network, layer, color):
    lines = []
    for node in network.nodes_where({'is_anchor': True}):
        start = network.node_attributes(node, 'xyz')
        residual = network.node_attributes(node, ['rx', 'ry', 'rz'])
        end = subtract_vectors(start, residual)
        lines.append({'start': start, 'end': end, 'arrow': 'end', 'color': color})
    compas_rhino.draw_lines(lines, layer=layer)


def draw_residuals(network, layer, color, tol):
    lines = []
    for node in network.nodes_where({'is_anchor': False}):
        start = network.node_attributes(node, 'xyz')
        residual = network.node_attributes(node, ['rx', 'ry', 'rz'])
        if length_vector(residual) < tol:
            continue
        end = add_vectors(start, residual)
        lines.append({'start': start, 'end': end, 'arrow': 'end', 'color': color})
    compas_rhino.draw_lines(lines, layer=layer)

Now we can use these function directly to implement our form-finding procedure to find equilibrium: As long as the residuals are below a threshold, no equilibrium is found yet, so we must iterate over and over again.

We limit the procedure to a maximum iteration number, so in case the equilibrium cannot be found, our script does not run forever.

In code this looks like this:

# ==============================================================================
# form finding
# ==============================================================================

# define maximum iterations and tolerance for residuals 
tol = 0.01
kmax = 100

# compute the residual forces in the original geometry
update_residuals(network)

for k in range(kmax):

    # stopping criteria
    R = network.nodes_attributes(['rx', 'ry', 'rz'], keys=free)
    res = sum(length_vector(r) for r in R)
    if res < tol:
        break
        
    # visualize dynamic process >
    # artist.draw_nodes(color={node: (255, 0, 0) for node in fixed})
    # artist.draw_edges()

    # draw_reactions(network, layer, (0, 255, 0))
    # draw_residuals(network, layer, (0, 255, 255), tol)

    # compas_rhino.rs.Redraw()
    # compas_rhino.wait()
    # < visualize dynamic process

    # update the geometry based on the previous residuals 
    update_geometry(network)
    # then recompute the residuals in the new geometry 
    update_residuals(network)

# ==============================================================================
# visualization
# ==============================================================================

artist.draw_nodes(color={node: (255, 0, 0) for node in network.nodes_where({'is_anchor': True})})
artist.draw_edges()

draw_reactions(network, layer, (0, 255, 0))
draw_residuals(network, layer, (0, 255, 255), tol)

In order to see the iteration steps, uncomment the lines in the loop:

e.x further exploration with non-uniform forces

Try what happens if you had a higher force in one of the edges:

network.add_edge(a, e, f=2)

This is what you get:

And this is how the iterations look like to get there:

It is not typically the straight path towards the equilibrium state as in our example with the uniform edge forces.

f. force densities

Instead of defining forces for the edges, we will now work with force densities. These are the ratio between the force and length of an edge. E.g. the same force causes a higher force density in a short edge than in a long edge, or the identical force density can be in one edge of high force and long length and one edge of little force and a short length. It can be considered as the density of woman/man-power pulling on a rope.

q:=FLq := \frac{F}{L}
# replace:
# network.update_dea(f=1)
# with:
network.update_dea(q=0.5)
            # replace:
            # F = network.edge_attribute(edge, 'f')
            # L = network.edge_length(*edge)
            # R[0] += F * (B[0] - A[0]) / L
            
            # with:
            Q = network.edge_attribute(edge, 'q')
            R[0] += Q * (B[0] - A[0])

We will see that the equation to compute the residual forces becomes much simpler and we must not divide by L.

If there is more than one free node, the numeric computation must be carried out with force densities to avoid a non-linear problem.

f.x further exploration with one longer edge

In order to understand how the force density vs. force input impacts the equilibrium in case of a different geometry, change the position of node b for both cases. 

# instead of:
# b = network.add_node(x=20, y=-10, z=20, is_anchor=True)
# try with a longer edge instead:
b = network.add_node(x=20, y=-10, z=20, is_anchor=True)

g. using a compas function fd_numpy

The procedure above is suitable for small data structures, however, as soon we deal with larger data structures it becomes highly inefficient. The same procedure in principle can be computed numerically in a more efficient manner.

The COMPAS framework already provides a numerical implementation of the force density method to compute the equilibrium of axial force networks with thefd_numpyfunction.

Using this function, this is how the hypar geometry of a cable net (now as Mesh data structure) looks before the form finding:

And this is how the hypar geometry of the cable net looks in its form found equilibrium state:

The upper equilibrium was found for uniform force densities of 1. Try increasing the force densities of the boundary edges to 5 to spread out the surface by 'pulling' on its boundaries:

for edge in mesh.edges_on_boundary():
    mesh.edge_attribute(edge, 'q', 5.0)

We will take a closer look at the logic of the script with our actual working geometry in step A2. Form Finding.

compas_fofin

The COMPAS package compas_fofin provides all functionalities for the form finding of cable-nets (also as flexible formworks). It enables their design with ease.

However, we don't want to just use a package, but rather understand the underlying logic. Thus we will implement it ourselves within the next three weeks!

PreviousA. Form FindingNextA1. Import

Last updated