I receive a MemError when I interpolate data with generator expressions and iterators from a regular grid on a mesh >2000 times

this is my first question here at stackoverflow, because I started scripting with Python3.

Application

I made a Python3 script that writes the load definition of a moveable heat source for a finite element simulation in LS-Dyna. As source I have a discretized 3D heat generation rate density (W/cm^3) field, the coordinates defining the finite element mesh and the position of the heat field center over time. As output I get a time dependent heating power sorted after the element number for each finite element. This works already for reasonable dimensions (200000 finite Elements, 3000 locations of the heat field, 400000 data points in the heat field).

Problem

For larger finite element meshes (4 000 000 Elements), I run out of memory (60GB RAM, python3 64Bit). To illustrate the problem further I prepared a minimal example which runs on its own. It generates some artificial test data, a finite element mesh how I use it (in reality it is not a regular grid) and an iterator for new locations for the heat application.

import numpy as np
import math
from scipy.interpolate import RegularGridInterpolator
def main():
    dataCoordinateAxes,dataArray = makeTestData()
    meshInformationArray = makeSampleMesh()
    coordinates = makeSampleCoordinates()
    interpolateOnMesh(dataCoordinateAxes,dataArray,meshInformationArray,coordinates)

def makeTestData():
    x = np.linspace(-0.02,0.02,300)
    y = np.linspace(-0.02,0.02,300)
    z = np.linspace(-0.005,0.005,4)
    data = f(*np.meshgrid(x,y,z,indexing='ij',sparse=True))
    return (x,y,z),data

def f(x,y,z):
    scaling = 1E18
    sigmaXY = 0.01
    muXY = 0
    sigmaZ = 0.5
    muZ = 0.005
    return weight(x,1E-4,muXY,sigmaXY)*weight(y,1E-4,muXY,sigmaXY)*weight(z,0.1,muZ,sigmaZ)*scaling

def weight(x,dx,mu,sigma):
    result = np.multiply(np.divide(np.exp(np.divide(np.square(np.subtract(x,mu)),(-2*sigma**2))),math.sqrt(2*math.pi*sigma**2.)),dx)
    return result

def makeSampleMesh():
    meshInformation = []
    for x in np.linspace(-0.3,0.3,450):
        for y in np.linspace(-0.3,0.3,450):
            for z in np.linspace(-0.005,0.005,5):
                meshInformation.append([x,y,z])
    return np.array(meshInformation)

def makeSampleCoordinates():
    x = np.linspace(-0.2,0.2,500)
    y = np.sqrt(np.subtract(0.2**2,np.square(x)))
    return (np.array([element[0],element[1],0])for element in zip(x,y))

The interpolation is then done in this function. I removed everything in the for loop to isolate the problem. In reality I export the load curve to a file in a specific format.

def interpolateOnMesh(dataCoordinateAxes,dataArray,meshInformationArray,coordinates):
    interpolationFunction = RegularGridInterpolator(dataCoordinateAxes, dataArray, bounds_error=False, fill_value=None)
    for finiteElementNumber, heatGenerationCurve in enumerate(iterateOverFiniteElements(meshInformationArray, coordinates, interpolationFunction)):
        pass
    return

def iterateOverFiniteElements(meshInformationArray, coordinates, interpolationFunction):
    meshDataIterator = (np.nditer(interpolationFunction(np.subtract(meshInformationArray,coordinateSystem))) for coordinateSystem in coordinates)
    for heatGenerationCurve in zip(*meshDataIterator):
        yield heatGenerationCurve

if __name__ == '__main__':
    main()

To identify the problem, I tracked the memory consumption over time. Memory Consumption over Time It seems the iteration over the result arrays consumes a considerable amount of memory.

Question

Is there a less memory consuming way to iterate over the datapoints without loosing too much performance? If not, I guess I will slice the mesh array in chunks and interpolate on those one by one.