C1.1 Dynamic equilibrium

a. Laws of motion

Newton's second law of motion states that the rate of change of momentum (mass . acceleration) of a body over time is proportional to the force applied on the body:

$$\textbf{F = m . a} \\~\\ v = \frac{dx}{dt} \\~\\ a = \frac{dv}{dt}$$

By further writing out the terms of this basic equation, we can describe the instantaneous displacement of each node of a structural network, under residual forces acting on it, per iteration step:

$$\begin{aligned} x_1 & = x_0 + dx\\ & = x_0 + (v_0 + dv) \cdot dt\\ & = x_0 + (v_0 + a \cdot dt) \cdot dt\\ & = x_0 + (v_0 + \frac{F}{m} \cdot dt) \cdot dt\\ \end{aligned}$$

Setting the following parameters:

$$\begin{aligned} & \text{initial velocity } v_0 = 0\\ & \text{mass } m = 1\\ & \text{time step } dt = 1\\ &\text{force } F = residual \end{aligned}$$

we can re-state the iterative nodal equilibrium equations from the above as:

$$\begin{aligned} x_1 & = x_0 + (0 + r \cdot 1) \cdot 1\\ & = x_0 + r \end{aligned}$$

This is again the simple updated residual equation that we wrote in our naive form finding method.

b. Damping

In C0.1 First Steps, the iterative procedure has been to update the XYZ coordinates of each free node with a factor (0.0 - 1.0) of the residuals at that node. This factor can be seen as a damping factor: as if the free nodes move through a viscous medium that exercises a force countering their displacements.

f_damping = 0.5

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 + f_damping * rx
        y1 = y0 + f_damping * ry
        z1 = z0 + f_damping * rz
        network.node_attributes(node, 'xyz', [x1, y1, z1])

The damping factor helps to avoid overshooting an equilibrium position and a potentialy infinite iteration of nodal coordinates oscilating around the equilibrium state.

DR: damped nodal oscillation around equilibrium position

The same concept is applied in dynamic relaxation: in order to avoid overshooting an equilibrium position, different damping strategies can be applied. Here, a constant viscous damping factor is applied to both the velocity and acceleration terms. The exact ratio between these two damping parameters can be optimized for the fastest convergence of the algorithm.

$$x = x_0 + (\textcolor{red}{c_1} \cdot v_0 + \textcolor{red}{c_2} \cdot a \cdot dt) \cdot dt$$

c. COMPAS implementation

The compas.numerical package has various dynamic relaxation functions available. We will be using the pure Python implementation, with no external linear algebra libraries (e.g. Numpy) for solvers.

from compas.numerical import dr

# ==============================================================================
# numerical data
# ==============================================================================

X = network.nodes_attributes('xyz')
R = network.nodes_attributes(['rx', 'ry', 'rz'])
P = network.nodes_attributes(['px', 'py', 'pz'])
Q = network.edges_attribute('q')


# ==============================================================================
# compute equilibrium
# ==============================================================================

X, Q, F, L, R = dr(X, edges, fixed, P, Q)

# update network
update_network()

d. Iteration control

In order to visualize the intermediate steps of the iterative algorithm, we can use one of the dr optional arguments which accepts a callable. The function provided here is called back at each iteration step. So if we provide a callback function which visualizes updated geometry and forces, each step shows the progress of the algorithm.

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

def callback_visualize(k, X, crits, args):
    if k % 2 == 0:
        # update nodal coordinates
        compas_rhino.wait()
        for node in network.nodes():
            index = node_index[node]
            network.node_attributes(node, 'xyz', X[index])

        # visualize updated geometry
        artist.draw_nodes(color={node: (255, 0, 0) for node in network.nodes_where({'is_anchor': True})})
        artist.draw_edges()
        draw_loads(network, layer, (255, 0, 0))
        compas_rhino.rs.Redraw()
        compas_rhino.wait()


# ==============================================================================
# compute equilibrium
# ==============================================================================

# calling dynamic relaxation function
X, Q, F, L, R = dr(X, edges, fixed, P, Q, callback=callback_visualize)