Is there a built-in function for interpolating a GridFunction/CoefficientFunction onto an n-dimensional grid for export to numpy, i.e. the inverse of VoxelCoefficient? There are a lot of nice built-in functions for evaluation of data, but there are limits, and I think at some point I will need to use numpy for analysis. I’ve searched tutorials and the forum and only found one user-written function to do this that I’m not completely sure about.

I ended up adapting a function posted by s.mokbel here for 3D with a scalar u, and a touch more robustness (np.nan for outside mesh, and uses maximum and minimum points in the mesh for the interpolation).

def create_np_data(gfu: GridFunction, mesh: Mesh, l: int, m: int, n: int):

“”"

Function to create numpy array from Gridfunction data.

Args:

gfu: Scalar gridfunction

mesh: The mesh used in this simulation.

l: Desired number of discrete cells (elements) in the uniform grid (numpy array)

in the x direction.

n: Desired number of discrete cells (elements) in the uniform grid (numpy array)

in the y direction.

m: Desired number of discrete cells (elements) in the uniform grid (numpy array)

in the z direction.

Returns:

Numpy array of field data U

“”"

# Obtain mesh dimensions from vertices.

vertices = np.array([vertex.point for vertex in mesh.vertices])

max_vertex = vertices.max(axis=0)

min_vertex = vertices.min(axis=0)

x0, y0, z0 = min_vertex

xn, yn, zn = max_vertex

# Initialize numpy arrays for field data.

output_U = np.zeros((l,m,n),dtype=float)

# Create uniform grid points based on mesh dimensions.

x_interp = np.linspace(x0,xn,l)

y_interp = np.linspace(y0,yn,m)

z_interp = np.linspace(z0,zn,n)

```
x_idx = 0
for p1 in x_interp:
y_idx = 0
for p2 in y_interp:
z_idx = 0
for p3 in z_interp:
try:
val_U = gfu(mesh(p1,p2,p3))
output_U[x_idx, y_idx, z_idx] = val_U
except Exception:
# Catch any points that might not be in the mesh,
# and set their values to np.nan
output_U[x_idx, y_idx, z_idx] = np.nan
z_idx += 1
y_idx += 1
x_idx += 1
return output_U
```