C0.3 Numerical Implementation

Objectives

We make the important transition from geometrical attributes to linear algebra / matrices (here represented as nested lists).

i. Attributes to lists

Imagine applying the equilibrium method we developed to a large structural network with many nodes and edges. The algorithm might still work, but the data storage and computational cost can quickly become problematic with the current naive implementation. That's where the numerical implementation comes in: by storing the data in nested lists (~ matrices), iteration over large amounts of data becomes manageable, and the form-finding methods can be applied efficiently. Firstly, some helper variables are introduced, and the data is implemented as (nested) lists:

n = network.number_of_nodes()

# mapping of node keys to contiguous node indices
node_index = {node: index for index, node in enumerate(network.nodes())}

# indices of fixed and free nodes
fixed[:] = [node_index[node] for node in fixed]
free[:] = [node_index[node] for node in free]

X = network.nodes_attributes('xyz')
R = network.nodes_attributes(['rx', 'ry', 'rz'])

# mapping of node index to indices of all its neighbors
i_nbrs = {node_index[node]: [node_index[nbr] for nbr in network.neighbors(node)] for node in network.nodes()}

# mapping of edge tuple to force densities
ij_force_densities = {}
for u, v in network.edges():
    i = node_index[u]
    j = node_index[v]
    force = network.edge_attribute((u, v), 'q')
    ij_force[i, j] = force

The numerical function inputs are lists that are ordered in the index and not the keys (vertex indices), so we must convert for the input the keys to indices and for the returns the indices to keys!

node_index = {node: index for index, node in enumerate(network.nodes())}

node_index generates a mapping of contiguous node indices

Secondly, the equilibrium helper functions are updated:

def update_network():
    for node in network.nodes():
        index = node_index[node]
        network.node_attributes(node, ['x', 'y', 'z'], X[index])
        network.node_attributes(node, ['rx', 'ry', 'rz'], R[index])

# replace:
# def update_residuals(network):
#     for node in network.nodes():
#         A = network.node_attributes(node, '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

#             Q = network.edge_attribute(edge, 'q')
#             R[0] += Q * (B[0] - A[0])
#             R[1] += Q * (B[1] - A[1])
#             R[2] += Q * (B[2] - A[2])
#         network.node_attributes(node, ['rx', 'ry', 'rz'], R)
# with:
def update_R():
    for i in range(n):
        R[i] = [0, 0, 0]
        a = X[i]
        for j in i_nbrs[i]:
            q = ij_force_densities[i, j]
            R[i][0] += q * (b[0] - a[0])
            R[i][1] += q * (b[1] - a[1])
            R[i][2] += q * (b[2] - a[2])

# and replace:
# 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])
# with:
def update_X():
    for i in range(n):
        if i in fixed:
            continue
        X[i][0] += 0.5 * R[i][0]
        X[i][1] += 0.5 * R[i][1]
        X[i][2] += 0.5 * R[i][2]

Finally, the solver routine is redefined with the new helper functions:

# ==============================================================================
# iterative equilibrium
# ==============================================================================

for k in range(kmax):
    if k % 10 == 0:
        if sum(length_vector(R[i]) for i in free) < tol:
            break
    
    # update nested lists of coordinates and residuals
    update_X()
    update_R()

The numerical implementation allows us to leverage efficient linear algebra solver algorithms in the backend. An explanation of linear algebra methods is out of scope of this workshop. It suffices to remember the importance of efficient and robust numerical implementation when scaling up the number of elements in a structural system.