I need to assign data from a mat file to a fem.Function, so that I can use it as the initial condition for my solver.
I don’t know how to ensure that the values from my mat file are read in in the correct order, as I understand that the ordering of vertices in the mesh doesn’t correspond to the ordering of the dofs in a function. So if I simply reshape my (7,7,3) data to be shape (147,), and assign these values to a function, the order of this reshaping is wrong. I know how to get the dofs of the FunctionSpace, but I’m not sure how to use this info.
I’ve looked at a lot of posts on the FEniCS legacy forum, in particular https://fenicsproject.org/qa/2715/coordinates-u_nodal_values-using-numerical-source-function/?show=2721#a2721 and https://fenicsproject.org/qa/3258/manually-setting-values-to-a-function/ , but they use legacy FEniCS parameters which don’t exist in FEniCSx.
I’d be very grateful for any pointers on this! Here is my MWE, I can provide the data mat file if it’s helpful.
import dolfinx
import numpy as np
import ufl
import scipy.io
from dolfinx import fem, mesh
from mpi4py import MPI
# Load in my data:
NTC1 = scipy.io.loadmat('testcase1/NTC1_testcase1.mat')
NTC1 = NTC1["NTC1"]
# Make mesh based on size of data:
L = NTC1.shape[0]-1
W = NTC1.shape[2]-1
mesh = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, L, W])],
[L, L, W], cell_type=mesh.CellType.hexahedron)
# Define function space
V = fem.FunctionSpace(mesh,("CG",1))
# Get the dofs of the FunctionSpace, I'm not actually using this currently:
local_range = V.dofmap.index_map.local_range
dofs = np.arange(*local_range)
# Reshaping the (7,7,3) data to be (147,), possibly the step where it all goes wrong
NTC1 = np.reshape(NTC1, (mesh.geometry.x.shape[0],));
# Define function (that we want to put NTC1 values into):
N1 = fem.Function(V)
N1.vector[:] = NTC1 # and this messes up the ordering
Thanks,
Rose
Hi @dokken, thanks for the quick answer! Unfortunately I’m quite new to FEniCS so I’m struggling to understand your code and apply it to my case. My main questions are:
-
How do I apply this to my case, ie to the case of having data in a matrix loaded in from a mat file, which I want to assign to the dofs of a fem.Function? I have tried the following, but the ordering of nodes is still getting mixed up. When I print out vertex_to_dof_map, it is just a list of the indices in increasing number (0,1,2,…,last_vertex) so it makes sense that my line N1.vector[vertex_to_dof_map[i]] = NTC1[i] does nothing.
-
Looking at your code, the main bit I don’t understand is the use of V0. What does V0, V0_to_V = V.sub(0).collapse() do – am I right in assuming it “collapses” the MixedElement space V onto the subspace of linear Lagrange elements V0? In that case, is it only necessary in the case of a MixedElement subspace (which I don’t require)? I also don’t understand the line bs = V0.dofmap.bs
Here is the code I mentioned in 1:
import dolfinx
from IPython import embed
import numpy as np
import pyvista
import ufl
import scipy.io
from dolfinx import fem, io, mesh, plot, nls, log
from dolfinx.io import gmshio
from mpi4py import MPI
from petsc4py import PETSc
# Load my data
NTC1 = scipy.io.loadmat('testcase1/NTC1_testcase1.mat')
NTC1 = NTC1["NTC1"]
# Make a mesh of corresponding dimensions
L = NTC1.shape[0]-1
W = NTC1.shape[2]-1
mesh = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, L, W])],
[L, L, W], cell_type=mesh.CellType.hexahedron)
# dokken's code on our FunctionSpace
el = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(mesh, el)
u = dolfinx.fem.Function(V)
V0, V0_to_V = V.sub(0).collapse()
dof_layout = V0.dofmap.dof_layout
num_vertices = mesh.topology.index_map(
0).size_local + mesh.topology.index_map(0).num_ghosts
vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
num_cells = mesh.topology.index_map(
mesh.topology.dim).size_local + mesh.topology.index_map(
mesh.topology.dim).num_ghosts
c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
for cell in range(num_cells):
vertices = c_to_v.links(cell)
dofs = V0.dofmap.cell_dofs(cell)
for i, vertex in enumerate(vertices):
vertex_to_dof_map[vertex] = dofs[dof_layout.entity_dofs(0, i)]
geometry_indices = dolfinx.cpp.mesh.entities_to_geometry(
mesh, 0, np.arange(num_vertices, dtype=np.int32), False)
x = mesh.geometry.x
bs = V0.dofmap.bs
for vertex, geom_index in enumerate(geometry_indices):
dof = vertex_to_dof_map[vertex]
for b in range(bs):
parent_dof = V0_to_V[dof*bs+b]
u.x.array[parent_dof] = x[geom_index, b]
u0 = u.sub(0).collapse()
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u0]) as vtx:
vtx.write(0.0)
# Applying the vertex_to_dof_map to put N values into f function:
NTC1 = np.reshape(NTC1, (mesh.geometry.x.shape[0],)); #reshape from matrix to vector of num_vertices length (=num_dofs length)
N1 = fem.Function(V)
for i in np.ndindex(vertex_to_dof_map.shape):
N1.vector[vertex_to_dof_map[i]] = NTC1[i]
Thanks,
Rose
Consider something like the following code:
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import fem, mesh, io
L = 2
W= 3
nx = 3
ny = 5
nz = 7
msh = mesh.create_box(comm=MPI.COMM_WORLD,
points=((0.0, 0.0, 0.), (L, L, W)), n=(nx-1, ny-1, nz-1),
cell_type=mesh.CellType.hexahedron,)
dx, dy, dz = L/(nx-1), L/(ny-1), W/(nz-1)
values = np.arange(nx*ny*nz).reshape(nx, ny, nz)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V = fem.FunctionSpace(msh, P1)
dof_coordinates = V.tabulate_dof_coordinates()
# dof index map to i,j,k
map_vec = np.zeros((len(dof_coordinates), 3), dtype=np.int32)
eps = 2e-14
for dof_index, coordinate in enumerate(dof_coordinates):
map_vec[dof_index] = [(coordinate[0]+eps)//dx, (coordinate[1]+eps)//dy, (coordinate[2]+eps)//dz]
u = fem.Function(V)
for i, map_i in enumerate(map_vec):
u.x.array[i] = values[map_i[0], map_i[1], map_i[2]]
with io.XDMFFile(msh.comm, "u.xdmf", "w") as xdmf:
xdmf.write_mesh(msh)
xdmf.write_function(u)
1 Like
Thank you so much, I really appreciate it!