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.
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.
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