Reputation: 49
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,
model.nodes is a dictionary containing all of the Node objects, keyed by their id
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
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,
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
Upvotes: 1
Views: 94
Reputation: 53089
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)
Upvotes: 2