Assembly and np._ix

import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

6.1. Assembly and np._ix#

Large matrices for PDE simulations are usually constructed from smaller matrices. Particularly in the finite element method, this is usually done element-by-element. The mesh is defined with two arrays, one containing the coordinates of the nodes, and another containing information on how the nodes are connected in elements. In a real finite element code, the element matrix would depend on thecoordinates, but below we work with hard-coded element matrix, so we do not need to define node coordinates. The purpose is just to illustrate how the global matrix is assembled on a system that is small enough to be able to print the whole thing.

Example code for finite element matrix assembly#

import numpy as np

# Connectivity matrix
elements = np.array([ 
    [0, 1, 4], 
    [1, 2, 4], 
    [2, 3, 4], 
    [3, 0, 4] 
    ])

# Element matrix (simplified)
K_elem = np.array([
        [2.0, -1.0, -1.0],
        [-1.0, 2.0, -1.0],
        [-1.0, -1.0, 2.0]
    ])

n_nodes = np.max(elements)+1

# METHOD 1: Naive approach with double loop over the element matrix
K_naive = np.zeros((n_nodes, n_nodes))
for elem_nodes in elements:
    for i in range(3):
        for j in range(3):
            K_naive[elem_nodes[i], elem_nodes[j]] += K_elem[i, j]

# METHOD 2: Efficient approach using numpy.ix_
K_efficient = np.zeros((n_nodes, n_nodes))
for elem_nodes in elements:
    K_efficient[np.ix_(elem_nodes, elem_nodes)] += K_elem    


print("Global matrix from naive implementation")
print(K_naive)

print("\nGlobal matrix from efficient implementation")
print(K_efficient)

As you can see when running the cell above, the two assembly methods lead to the same matrix.

Questions

The global matrix contains four 0’s. Can you explain their position?

Solution

The entries on the global matrix that are related to nodes that are not connected in any element remains zero. In the simple mesh, that is the case for nodes 0 and 2, and for nodes 1 and 3. The zero entries are therefore at the positions K[0,2], K[2,0], K[1,3] and K[3,1].

The diagonal of the global matrix contains all 4’s and one 8. Can you explain why the values are as they are?

Solution

The element matrix is every time added to the global matrix. The element stiffness has 2’s on the diagonal, so the value on the diagonal in the global matrix will be equal to 2 times the number of elements that contain the corresponding nodes. Nodes 0, 1, 2 and 3 are in two elements each, while node 4 is in four elements.

Vector assembly#

In the context of this lesson it is good to recall that for indexing in vectors, we do not need np._ix. As an example, consider that we are just interested in the diagonal of the global matrix. This can be represented as a vector, which we can assemble by every time adding the diagonal vector of the element matrix to the correct entries.

d_elem = np.diag(K_elem)
d_global = np.zeros(n_nodes)

for elem_nodes in elements:
    d_global[elem_nodes] += d_elem

print("Global vector assembled without loop over nodes and without np._ix")
print(d_global)

Attribution

This chapter is written by Frans van der Meer. Find out more here.

Contents