Networks of bars and nodes

Posted on July 23rd, 2014, by tom in Uncategorized

Many calculations in structural design involve networks of bars (or branches, edges, …) and nodes (or vertices, …). Essential to these calculations is information about how the elements of the network are connected. This connectivity can be described with a matrix.

Connectivity matrix

Each row in a connectivity matrix corresponds to a bar of the network and contains zeros in all columns, except in the columns of the nodes to which it is connected. There the entries are ‘-1’ and ‘+1’ indicating the direction of the edge:

\begin{equation}
C_{ij} = \left\{
\begin{aligned}
+1 & \text{ if vertex } i \text{ is the head of edge } j\\
-1 & \text{ if vertex } i \text{ is the tail of edge } j\\
0 & \text{ otherwise}
\end{aligned}
\right.
\end{equation}

The direction of the edges can be chosen arbitrarily. A typical convention is to direct all edges from the vertex with the lower index to the one with the higher index. Constructing a connectivity matrix in Python is easy if you have list of edge vertex pairs.

import numpy as np

edges = [[0, 4], [1, 4], [2, 4], [3, 4]]
vertices = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0.5, 0.5, 0]]

m = len(edges)
n = len(vertices)

C = np.zeros((m, n))
for i in range(m):
    u, v = edges[i]
    C[i, u] = -1
    C[i, v] = +1

# print C
array([[-1.,  0.,  0.,  0.,  1.],
      [ 0., -1.,  0.,  0.,  1.],
      [ 0.,  0., -1.,  0.,  1.],
      [ 0.,  0.,  0., -1.,  1.]])

Once you have the connectivity matrix, it is easy to compute some of the properties of the elements of the network.

xyz = np.array(vertices)

# coordinate difference vectors
uvw = C.dot(xyz)

# edge lengths
l = np.sqrt(np.sum(np.square(uvw), axis=1))
# or, since operations like ** are elementwise
l = np.sum(uvw**2, axis=1)**0.5

# make l into a column vector
l = l.reshape((-1, 1))

# normalized coordinate difference vectors
unit = np.divide(uvw, np.tile(l, 3))
# or, since numpy does something called 'broadcasting'
unit = np.divide(uvw, l)
# or, since operations like / are elementwise
unit = uvw / l

Laplacian matrix, Degree matrix, Adjacency matrix

The laplacian matrix of a network can be derived from its connectivity matrix.

L = C.T.dot(C)

# print L
array([[ 1.,  0.,  0.,  0., -1.],
      [ 0.,  1.,  0.,  0., -1.],
      [ 0.,  0.,  1.,  0., -1.],
      [ 0.,  0.,  0.,  1., -1.],
      [-1., -1., -1., -1.,  4.]])

The laplacian matrix contains on its diagonal the degree of each vertex of the network.

degree = L.diagonal()

# print degree
array([ 1.,  1.,  1.,  1.,  4.])

The laplacian matrix is related to the degree matrix and the adjacency matrix in the following way:

\begin{equation}
\mathbf{L} = \mathbf{D} – \mathbf{A}
\end{equation}

Therefore, the adjacency matrix can be derived from the laplacian matrix like this

D = np.diagflat(degree)
A = D - L

# print A
array([[ 0.,  0.,  0.,  0.,  1.],
      [ 0.,  0.,  0.,  0.,  1.],
      [ 0.,  0.,  0.,  0.,  1.],
      [ 0.,  0.,  0.,  0.,  1.],
      [ 1.,  1.,  1.,  1.,  0.]])

Sparse matrices

Text will follow later…

from scipy.sparse import coo_matrix, diags

m = len(edges)
data = np.array([-1]*m + [1]*m)
rows = np.array(range(m) + range(m))
cols = np.array([edge[0] for edge in edges] + [edge[1] for edge in edges])
C = coo_matrix((data, (rows, cols))).tocsr()
L = C.transpose().dot(C)
D = diags([L.diagonal()], [0])
A = D - L

BRG blog