Convert a FEM solution function into an array

dolfinx
,
1

Hi, we would like to input our Fenicsx FEM solution into a NN.

Here is the domain and mesh parameters :

length, height = 1.0, 1.0   # Physical domain: 1 x 1
Nx, Ny = 64, 64            # Mesh resolution
domain = create_rectangle(
    MPI.COMM_WORLD,
    [np.array([0.0, 0.0]), np.array([length, height])],
    [Nx, Ny],
    cell_type=CellType.quadrilateral,
)

dim = domain.topology.dim

The solution is a (2D) displacement field : u_sol obtained as follows :

problem = fem.petsc.LinearProblem(
    a, L, u=u_sol, bcs=bcs,
    petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
problem.solve()

For that we would need to obtain an array with the following shape (2,64,64), one 2d displacement array per direction.
In order to do so we interpolated u_sol doing element wise interpolation (see code bellow):

# Create a DG-0 function space for interpolation (element-wise interpolation)
V_DG = fem.functionspace(domain, ("DG", 0, (dim,)))
# Create a new function in the DG space
u_dg_interpolated = fem.Function(V_DG, name="Interpolated Solution")
# Project the solution onto the DG space
u_dg_interpolated.interpolate(u_sol)
elem_values=u_dg_interpolated.x.array[:]
elem_val_x1 =elem_values[::2].reshape(64,64)
elem_val_y1 = elem_values[1::2].reshape(64,64)

However, when plotting these arrays it does not match the FEM solution visualized on paraview. We have also reshaped the array in other ways assuming the displacement are arranged differently (first 4096 values corresponds to x direction, last 4096 to y direction) but none of that seems to work. We would really appreciate some help to obtain an array that correctly represent the displacements.

Thank you in advance !

2

See Reordering nodal values in accordance with meshgrid created using numpy - #5 by AIJAZ_NAZIR