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
