Set values of Expression using a numpy.ndarray


106
views
1
4 months ago by
I want to solve a differential equation where I need to compute its source term "outside of FEniCS functionality". Assume we have computed the source term and stored its values in a numpy.ndarray. It is an array with 3 indices for the 3 spatial variables where the first index corresponds to the x coordinate, the second to the y-coordinate, and the third to the z-coordinate. I now need to reorganize this array to be 1d where the order of the elements corresponds to the mesh coordinates given by FEniCS. Of course I could write a function that maps each (x,y,z) value to an index in the mesh coordinates. I want to know if there is a better option. Heres my (not correctly working) attempt:
from fenics import *
import numpy as np

nx = ny = nz = 20
xmin = ymin = zmin = -1
xmax = ymax = zmax = 1
mesh = BoxMesh(Point(xmin, ymin, zmin), Point(xmax, ymax, zmax), nx, ny, nz)
coords = mesh.coordinates()
V = FunctionSpace(mesh, 'P', 1)

#########################################
#########################################
#########################################

# Define source term with an fenics.Expression which is what I like to emulate
I = Expression('1.0/(1+pow(x[2], 2))*exp(-pow(x[0], 2) - pow(x[1], 2))', degree=2)
I_n = interpolate(I, V)

# Define source term from a non-FEniCS computation which is what I need to do
x = np.linspace(xmin, xmax, nx+1)
y = np.linspace(ymin, ymax, ny+1)
z = np.linspace(zmin, zmax, nz+1)
I2 = np.zeros((nx+1, ny+1, nz+1))
for i in range(nx+1):
    for j in range(ny+1):
        for k in range(nz+1):
            I2[i, j, k] = 1/(1+z[k]**2)*np.exp(-x[i]**2-y[j]**2)

# Reshape I2 so that its components correspond to coordinates given by coords
I2 = I2.reshape((nx+1)*(ny+1)*(nz+1)) # this is most certainly wrong

# Set array of FEniCS expression to the numpy array
I_n2 = Function(V)
I_n2.vector().set_local(I2)
as_backend_type(I_n.vector()).vec().ghostUpdate()
print(np.max(np.abs(I_n2.vector().get_local()-I_n.vector().get_local())))

#########################################
#########################################
#########################################​

Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »