C0.1 First Steps

Iterative equilibrium

Objectives

We will gain a basic understanding of equilibrium methods by building a naive form finding algorithm step-by-step. As always, we compute all data steps with COMPAS, but for these steps we will visualize the resulting geometry inside the Rhino3D environment.

a. Network

First, we will create a network data structure with a geometry that is far from its equilibrium state. In this simple case with 5 nodes, we can generate the network from scratch.

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

  • is_anchor: False - to indicate that a node is not anchored (i.e. free in XYZ).

  • rx: 0.0, ry: 0.0, rz: 0.0 - no nodal residual nodal forces (yet).

  • f: 1 - uniform, unit edge force.

• The network will consist of five nodes (identified as a to e) of which the leaves are anchored, and of edges that all connect to one 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 = "DF21_C0::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()

We will place all the geometry in Rhino for this section C0. Principles inside the parent layer DF21_C0 and all the subsequent steps in sublayers. With :: you can indicate the sublayer.

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 here 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 vector sum of all forces in the edges (and any external loads) neighboring a node is the residual forces in that node; these are basically the leftover force in the free nodes (shown in cyan). For a static structure to be in equilibrium, there can be no leftover forces - except at the anchors (i.e. supports): here the residual forces are the reaction forces (shown in green).

So if all the edge forces acting on the free node were in equilibrium, the sum of these forces would be equal to zero in all x-, y- and z-directions:

node with edge forces in equilibrium

We can write out the equilibrium in each independent coordinate direction:

$$F_{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 randomly generated geometry as this one, it is very unlikely that all forces cancel out, and the edge forces would sum up to a residual force R (in red):

node with edge forces with residual force (in red)

$$F_{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 way through summation:

$$\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. For each of the coordinates, the direction cosine is given by the difference in coordinates divided by the edge length:

$$\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 first 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 we must separate the forces into its x-, y-, and z-components. We store the sum of residuals to the nodes by updating their 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)

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 their orientation flipped:

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

lines = []
for node in network.nodes():
    start = network.node_attributes(node, 'xyz')
    residual = network.node_attributes(node, ['rx', 'ry', 'rz'])
    if network.node_attribute(node, 'is_anchor'):
        end = subtract_vectors(start, residual)
        color = (0, 255, 0)
    else:
        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 given threshold), we want the node to 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 overshoot the equilibrium coordinates, we will apply a factor of 0.5 on the residuals for the displacement of the node.

Remember to recompute the new residual force after updating the geometry as it will have changed! Since we are re-using the code to calculate the residuals, we can pack it into a function.

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


# ==============================================================================
# move free nodes in direction of residual
# ==============================================================================
# compute residuals in the initial geometry
update_residuals(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])


# ==============================================================================
# recompute residuals in the updated geometry
# ==============================================================================

update_residuals(network)

e. Equilibrium

We must repeat the previous steps until equilibrium is found, meaning the residual forces in each free node are zero (or below a given threshold).

In order to keep our script clean and not to repeat code, we pack as much as possible the repeatable code blocks into 'helper' 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 functions directly to implement our form-finding procedure to find equilibrium. as long as the residuals are below a certain threshold, no equilibrium is found yet, so we must iterate over and over again.

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

Solving Zeno's paradox, or "close enough"

In code this looks as follows:

# ==============================================================================
# iterative equilibrium > NEW
# ==============================================================================
# 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):
    R = network.nodes_attributes(['rx', 'ry', 'rz'], keys=free)
    res = sum(length_vector(r) for r in R)
    # stopping criterion
    if res < tol:
        break

    # update the geometry based on the residuals of the previous step
    update_geometry(network)
    # 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)

f. Iteration control

In order to see the intermediate iteration steps, the visualisation can be moved inside the body of the for loop. In order to reduce the amount of visual information output, the functions can be called only every multiple of a number of iteration steps. The same can be done for checking the stopping criterion.

# ==============================================================================
# iterative equilibrium
# ==============================================================================
# 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):
    R = network.nodes_attributes(['rx', 'ry', 'rz'], keys=free)
    res = sum(length_vector(r) for r in R)
    # stopping criterion
    if k % 10 == 0:
        if res < tol:
            break

    # visualize the iterative process
    if k % 2 == 0:
        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()

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

g. External loads

External loads applied to the nodes can be treated similarly to edge forces. In principle, external loads may also be updated at each iteration step for the new geometry. For example, think of loads that change their amplitude based on the surrounding surface area of the node; or pneumatic loads which are always oriented along the local normals of a surface. In this simple example, the loads are held constant for all of the steps.

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

def update_residuals(network):
    for node in network.nodes():
        A = network.node_attributes(node, 'xyz')
        
        # external loads P as contribution to nodal residuals
        # R = P + F * cos(XYZ)
        R = network.node_attributes(node, ['px', 'py', 'pz'])

        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 draw_loads(network, layer, color):
    lines = []
    for node in network.nodes():
        start = network.node_attributes(node, 'xyz')
        load = network.node_attributes(node, ['px', 'py', 'pz'])
        end = add_vectors(start, load)
        lines.append(
            {'start': start,
             'end': end,
             'arrow': 'end',
             'color': color})
    compas_rhino.draw_lines(lines, layer=layer)


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

# add external load vector to the free node
network.node_attributes(e, ['px', 'py', 'pz'], [0, 0, -3])


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

compas_rhino.clear()
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)
draw_loads(network, layer, (255, 0, 0))

Try it yourself

• Taking the script C01_e_equilibirum, what would be the easiest way to ensure that the free node lies inside the plane defined by nodes a, b, and c? (hint: 4 independent forces can fully equilibrate a node in 3D space, 3 independent forces can equilibrate a node in 2D).