Importing and interpolating an initial condition from a .mat file

I/O
1

Good day everybody,

I am a beginner at using FEniCS, I use the latets version of FEniCS installed with Docker on Windows 10. I have the following problem:
I would like to solve a transient diffusion PDE with homogeneous Dirichlet and Neumann BCs with an arbitrary initial condition that I import from a .mat file. In the end, the goal is to export the solution back to Matlab for postprocessing. At the end is a (hopefully) minimal example of my code.

My question is about the initial condition:
Based on an older answer to a similar question (From an older question (https://fenicsproject.org/qa/465/interpolate-function-given-matrix-on-finite-element-basis/) I tried to define my initial condition in the same way, but the resulting initial state was still wrong, looking scrambled. The problem seems to be the mapping of the content of the .mat file to correct vertices of the mesh, but i am unaware how to properly do that in the current version, since the solutions in the linked question (parameters['reorder_dofs_serial'] = False and psi.vector()[:] = phi_vector[V.dofmap().vertex_to_dof_map(mesh)]) both didn’t work for me, the end result after exporting the data still looked scrambled.

In the example below I pretty much only set the IC and then export it, here I would like to know how to properly set the IC so that the exported file looks the same as the imported one.

Working example:

from __future__ import print_function
from fenics import *
import numpy as np
from scipy.io import loadmat, savemat

T = 0.5            # final time
X = 0.5
Y = 1.0

# Subdomains
class Inlet(SubDomain):
     Y = Y;
     def inside(self, x, on_boundary):
         return x[1] > Y - DOLFIN_EPS and on_boundary

class Outlet(SubDomain):
     def inside(self, x, on_boundary):
         return x[1] < DOLFIN_EPS and on_boundary

class Walls(SubDomain):
     def inside(self, x, on_boundary):
         return on_boundary

# Create mesh and define function space
nx = ny = 5
mesh = RectangleMesh(Point(0,0),Point(X,Y),nx,ny)
V = FunctionSpace(mesh, 'CG', 1)

# Subdomains
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1);

walls = Walls()
inlet = Inlet()
outlet = Outlet()

u_Diri = Constant('0')
bc_Inlet = DirichletBC(V, u_Diri, inlet)
bc_Outlet = DirichletBC(V, u_Diri, outlet)

bcs = [bc_Inlet,bc_Outlet]

# Define initial value

u_Ini = Function(V)

# Load Initial State from .mat
ini_cont = loadmat('example.mat')
ini = ini_cont['M']
ini_vector = (ini[:,:]).reshape((nx+1)*(ny+1))
u_Ini.vector()[:] = ini_vector[vertex_to_dof_map(V)]

u_n = interpolate(u_Ini, V)

t = 0

hdf = HDF5File(mesh.mpi_comm(),'InitialState.h5','w')
# Export the first step
hdf.write(u_n,'fun',t)

# Calculations are excluded here...
...
...
...
# Export the coordinates in a separate .csv file
Coord = V.tabulate_dof_coordinates()
Coord.resize((V.dim(), 2))
x = Coord[:,0]
z = Coord[:,1]
np.savetxt("Results/Coords.csv", np.c_[x, z], delimiter=",")