Do I really need to re-read your code and compare it to the one I created above?
Why not simply use the code I supplied, i.e. the function:
There is a reason for me adding all the steps above, hich is that you need to go through:
Of course, if you have a higher order mesh, and want to map nodes from say a midpoint of a facet to the dof at a midpoint of a facet, something slightly different should be done.
"""
Compute local dof indices of global input node indices
Author: Jørgen S. Dokken and Henrik N.T. Finsberg
SPDX-License-Identifier: MIT
"""
from mpi4py import MPI
import dolfinx
import numpy as np
def input_indices_to_dof(V, global_dofs):
mesh = V.mesh
mesh.topology.create_connectivity(0, mesh.topology.dim)
# Create a mapping from global to local node indices
num_nodes_global = mesh.geometry.index_map().size_global
local_to_input_global = mesh.geometry.input_global_indices
global_to_local = np.full(num_nodes_global, -1, dtype=np.int32)
global_to_local[local_to_input_global] = np.arange(len(local_to_input_global), dtype=np.int32)
geom_layout = mesh.geometry.cmap.create_dof_layout()
V_layout = V.dofmap.dof_layout
for i in range(mesh.topology.dim):
num_entities_per_cell = dolfinx.cpp.mesh.cell_num_entities(
mesh.topology.cell_type, i
)
for j in range(num_entities_per_cell):
gl = geom_layout.entity_dofs(i, j)
vl = V_layout.entity_dofs(i, j)
assert vl == gl, "Mismatch in dof layout"
geom_to_dof_layout = np.empty(V.dofmap.index_map.size_local+ V.dofmap.index_map.num_ghosts, dtype=np.int32)
geom_to_dof_layout[mesh.geometry.dofmap.reshape(-1)] = V.dofmap.list.reshape(-1)
# Remove all nodes that are not on process
local_nodes = global_to_local[global_dofs]
local_nodes = local_nodes[local_nodes >= 0]
return geom_to_dof_layout[local_nodes]
if __name__ == "__main__":
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
# Select a random set of input nodes for illustration
np.random.seed(0)
num_nodes_global = mesh.geometry.index_map().size_global
input_nodes = np.random.randint(0, num_nodes_global, 10)
print(MPI.COMM_WORLD.rank, input_nodes)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
dofs = input_indices_to_dof(V, input_nodes)
x_coord = V.tabulate_dof_coordinates()
# Sanity check
num_nodes_global = mesh.geometry.index_map().size_global
local_to_input_global = mesh.geometry.input_global_indices
global_to_local = np.full(num_nodes_global, -1, dtype=np.int32)
global_to_local[local_to_input_global] = np.arange(len(local_to_input_global), dtype=np.int32)
local_geom = global_to_local[input_nodes]
local_geom = local_geom[local_geom >= 0]
for local_node, dof in zip(local_geom, dofs):
assert np.allclose(mesh.geometry.x[local_node], x_coord[dof])
print(mesh.geometry.x[local_node,], x_coord[dof])
This will only work if you use the same finite element to describe the geometry of your mesh as you would use for the function space you want to store data in.