I don’t have time to read that user manual.
Assuming that what is in NSET describes the original nodes in the mesh you would like to constrain, you need to:
- Create a map from the input node index to a vertex on your process, this is simply a matter of inverting
mesh.geometry.input_global_indices. - Once you have a map from the global indices (in your case a map from
97, 96, 95, ...to what is local on the process), you remove all the results that are -1. - You then map from geometry to topology indices, and then from topology to dof indices.
I’ve illustrated the procedure below using a set of random input nodes (based on @Henrik_Nicolay_Finsb vertorized vertex_to_dofmap:
"""
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
# Compute map from vertex to dof
def vertex_to_dof_map_vectorized(V):
mesh = V.mesh
num_vertices_per_cell = dolfinx.cpp.mesh.cell_num_entities(
mesh.topology.cell_type, 0
)
dof_layout2 = np.empty((num_vertices_per_cell,), dtype=np.int32)
for i in range(num_vertices_per_cell):
var = V.dofmap.dof_layout.entity_dofs(0, i)
assert len(var) == 1
dof_layout2[i] = var[0]
num_vertices = (
mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
)
c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
assert (
c_to_v.offsets[1:] - c_to_v.offsets[:-1] == c_to_v.offsets[1]
).all(), "Single cell type supported"
vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
vertex_to_dof_map[c_to_v.array] = V.dofmap.list[:, dof_layout2].reshape(-1)
return vertex_to_dof_map
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)
# Compute map from geometry node to vertex
num_vertices_local = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
vertices_to_geometry = dolfinx.mesh.entities_to_geometry(mesh, 0, np.arange(num_vertices_local, dtype=np.int32)).reshape(-1)
num_nodes_local = mesh.geometry.index_map().size_local + mesh.geometry.index_map().num_ghosts
geometry_to_vertices = np.full(num_nodes_local, -1, dtype=np.int32)
geometry_to_vertices[vertices_to_geometry] = np.arange(len(vertices_to_geometry), dtype=np.int32)
# Compute map from vertex to dof
v_to_d = vertex_to_dof_map_vectorized(V)
# Remove all nodes that are not on process
local_nodes = global_to_local[global_dofs]
local_nodes = local_nodes[local_nodes >= 0]
local_vertices = geometry_to_vertices[local_nodes]
local_dofs = v_to_d[local_vertices]
return local_dofs
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])