Search code examples
pythonnumpydictionaryfinite-element-analysis

Looking for a better algorithm or data structure to improve conversion of connectivity from ID's to indices


I'm working with Python 3.6.2 and numpy.

I'm writing code to visualize a finite element model and results.

The visualization code requires the finite element mesh nodes and elements to be identified by indices (starting a zero, no gaps) but the input models are based on ID's and can have very large gaps in the ID space.

So I'm processing all of the nodes and elements and changing them to use indices instead of ID's.

The nodes are

First step is to process the array of nodes and node coordinates. This comes to me sorted so I don't specifically have to do anything with the coordinates - I just use the indices of the nodal coordinate array. But I do need to then redefine the connectivity of the elements to be index base instead of ID based.

To do this, I create a dictionary by iterating over the array of node ids and adding each node to the dictionary using it's ID as the key and its index as the value

In the following code fragment,

  1. model.nodes is a dictionary containing all of the Node objects, keyed by their id

  2. nodeCoords is a pre-allocated numpy array where I store the nodal coordinates for later use in visualization. It's the indices of this array that I need to use later to redefine my elements

  3. nodeIdIndexMap is a dictionary that I populate using the Node ID as the key and the index of nodeCoords as the value

Code:

nodeindex=0
node_id_index_map={}
for nid, node in sorted(model.nodes.items()): 
    nodeCoords[nodeIndex] = node.xyz   
    nodeIdIndexMap[nid] = nodeIndex
    nodeIndex+=1 

Then I iterate over all of the elements, looking up each element node ID in the dictionary, getting the index and replacing the ID with the index.

In the following code fragment,

  1. tet4Elements is a dictionary containing all elements of type tet4, keyed using the element id
  2. n1, n2, n3 and n4 are pre-allocated numpy arrays that hold the element nodes
  3. element.nodes[n].nid gets the element node ID
  4. n1[tet4Index] = nodeIdIndexMap[element.nodes[0].nid looks up the element node ID in the dictionary created in the previous fragment, returns the corresponding index and stores it in the numpy array

Code:

tet4Index = 0
for eid, element in tet4Elements.items():
    id[tet4Index] = eid
    n1[tet4Index] = nodeIdIndexMap[element.nodes[0].nid]
    n2[tet4Index] = nodeIdIndexMap[element.nodes[1].nid]
    n3[tet4Index] = nodeIdIndexMap[element.nodes[2].nid]
    n4[tet4Index] = nodeIdIndexMap[element.nodes[3].nid]
    tet4Index+=1 

The above works, but it's slow......It takes about 16 seconds to process 6,500,000 tet4 elements (each tet4 element has four nodes, each node ID has to be looked up in the dictionary, so that's 26 million dictionary lookups in a dictionary with 1,600,000 entries.

So the question is how to do this faster? At some point I'll move to C++ but for now I'm looking to improve performance in Python.

I'll be grateful for any ideas to improve performance.

Thanks,

Doug


Solution

  • With the numbers you are quoting and reasonable hardware (8GB ram) the mapping can be done in less than a second. The bad news is that getting the data out of the original dicts of objects takes 60 x longer at least with the mock objects I created.

    # extract 29.2821946144104 map 0.4702422618865967
    

    But maybe you can find some way of bulk querying your nodes and tets?

    Code:

    import numpy as np
    from time import time
    
    def mock_data(nn, nt, idf):
        nid = np.cumsum(np.random.randint(1, 2*idf, (nn,)))
        nodes = np.random.random((nn, 3))
        import collections
        node = collections.namedtuple('node', 'nid xyz')
        tet4 = collections.namedtuple('tet4', 'nodes')
        nodes = dict(zip(nid, map(node, nid, nodes)))
        eid = np.cumsum(np.random.randint(1, 2*idf, (nt,)))
        tet4s = nid[np.random.randint(0, nn, (nt, 4))]
        tet4s = dict(zip(eid, map(tet4, map(lambda t: [nodes[ti] for ti in t], tet4s))))
        return nodes, tet4s
    
    def f_extract(nodes, tet4s, limit=15*10**7):
        nid = np.array(list(nodes.keys()))
        from operator import attrgetter
        ncoords = np.array(list(map(attrgetter('xyz'), nodes.values())))
        tid = np.array(list(tet4s.keys()))
        tnodes = np.array([[n.nid for n in v.nodes] for v in tet4s.values()])
        return nid, ncoords, tid, tnodes, limit
    
    def f_lookup(nid, ncoords, tid, tnodes, limit):
        nmx = nid.max()
        if nmx < limit:
            nlookup = np.empty((nmx+1,), dtype=np.uint32)
            nlookup[nid] = np.arange(len(nid), dtype=np.uint32)
            tnodes = nlookup[tnodes]
            del nlookup
        else:
            nidx = np.argsort(nid)
            nid = nid[nidx]
            ncoords = ncoords[nidx]
            tnodes = nid.searchsorted(tnodes)
        tmx = tid.max()
        if tmx < limit:
            tlookup = np.empty((tmx+1,), dtype=np.uint32)
            tlookup[tid] = np.arange(len(tid), dtype=np.uint32)
        else:
            tidx = np.argsort(tid)
            tid = tid[tidx]
            tnodes = tnodes[tidx]
        return nid, ncoords, tid, tnodes
    
    data = mock_data(1_600_000, 6_500_000, 16)
    t0 = time()
    data = f_extract(*data)
    t1 = time()
    f_lookup(*data)
    t2 = time()
    print('extract', t1-t0, 'map', t2-t1)