Convert FEniCS solution and coordinates to (numpy) 2D arrays

I/O
1

Hello,

I am trying to convert FEniCS 2D field (scalar field or vector field) solutions and the corresponding x and y coordinates as 2D arrays in to separate Numpy arrays, got a bit confused by the conversion. Following is an example I am practicing on:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2)
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds
u = Function(V)
solve(a == L, u, bc)

in this case, solution u is a 2D scalar field, how should I convert it into three 2D numpy arrays T, X, Y such that each contains the field values, x coordinates, y coordinate, respectively?

Also, for the case when u is a 2D vector field of temperature, how should I convert it into four 2D numpy arrays Ux, Uy, X, Y?

Thank you!

1 Like
2

for scalar fields, use

xyz = V.tabulate_dof_coordinates()
x = xyz[:,0]
y = xyz[:,1]
array = u.vector().get_local()

For vector function spaces use

V0 = V.sub(0).collapse()
V1= V.sub(1).collapse()
ux, uy = u.split(deepcopy=True)

and then the same strategy as for scalar fields on each component.

4 Likes
4

Hello dokken, the above method gives 1D arrays only, I am trying to turn them into grid-like 2D arrays

5

See: Proper coordinates for eval()-method

6

Hi, I came across exactly the same issue as yours, do you still remember how you solved it?

7

fenics or fenicsx

In fenicsx I did the following and it worked

solver.solve(w)


# saving data 
tol = 0.0001  # avoid hitting outside of domain
den = 20 # take less points because the loop below is slow
x = np.linspace(tol,1 - tol,den)
y = np.concatenate([ np.linspace(tol,ep,np.int_(np.rint(bl_den*den*ep))),
               np.linspace(ep + 1/den, 1-ep - 1/den ,np.int_(np.rint(den*(1-2*ep - 2/den)))),
               np.linspace(1-ep,1 - tol, np.int_(np.rint(bl_den*den*ep)))  ] )




X,Y = np.meshgrid(x,y)
coords = np.array([ X.ravel(),Y.ravel(), np.zeros(X.ravel().shape)])

from dolfinx import geometry

bb_tree = geometry.BoundingBoxTree(msh, msh.topology.dim)
u_values = []
p_values = []
r_values = []
cells = []
points = []
# find cells whose bounding box collide with the points
cell_candidates = geometry.compute_collisions(bb_tree, coords.T)
# chooseone of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(msh, cell_candidates, coords.T)

for i, coord in enumerate(coords.T):
    if len(colliding_cells.links(i)) >0:
        points.append(coord)
        cells.append(colliding_cells.links(i)[0])

points = np.array(points,dtype = np.float64)
velocity_values = w.sub(0).eval(points, cells)
pressure_values = w.sub(1).eval(points, cells)
radius_values= w.sub(2).eval(points, cells)

sol = np.hstack((velocity_values, pressure_values, radius_values))



mesh_name = str( (den ,y.shape[0]) ) + '.dat'

# removing z = 0 column
coords = np.array([ X.ravel(),Y.ravel()]) 
np.savetxt(mesh_name, coords)
np.savetxt('sol.dat', sol)

While in legacy fenics

# Solve problem

solve(F == 0, w, bcs)

# Plot solutions
(u, p,r) = w.split()

xy = np.array(mesh.coordinates())


x = xy[:,0]
y = xy[:,1]

uw = np.array([u(pt) for pt in xy])
P_ = np.array([p(pt) for pt in xy])
R_ = np.array([r(pt) for pt in xy])


X = x.reshape(161,101)
Y = y.reshape(161,101)
u_ = uw[:,0]
w_ = uw[:,1]

U = u_.reshape(161,101 )
W = w_.reshape(161,101 )
P = P_.reshape(161,101 )
R = R_.reshape(161,101 )

where I reshaped the arrays based on the size of my mesh for post-processing in the legacy fenics example. Hope this helps.

8

Hi, thank you for your help! I’m using fenics, but when I ran w.split(), there are no return values and cannot be unpacked into three variables

9

Is the solver working? How is w defined.
I used split because w in my code is a mixed function space.

V = VectorElement("Lagrange", triangle, 2) # [ u,w] velocities
Q = FiniteElement("Lagrange", triangle, 1) # pressure
S = FiniteElement("Lagrange", triangle, 2) # free boundary problem

W = FunctionSpace(mesh, MixedElement([V, Q, S]))

# Define Function and TestFunction(s)
w = Function(W)

It depends on how you’ve defined w.

10

Thank you so much! I’ve solved the problem!

Furthermore, have you ever tried a Navier-Stokes problem in Fenics. I’ve got outputs of a 2D velocity field and V = FunctionSpace(mesh, “Lagrange”, 2). In this situation, I cannot extract the data into numpy array.

11

Should V be a vector function space
If it’s a funtionsapce you don’t need to use split