Extract the coordinates of the boundary points and the value of a function on them

dolfinx
4

Consider the following:

from mpi4py import MPI
import dolfinx
import numpy as np
import ufl
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 16, 16)
v_cg1 = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(mesh, v_cg1)
f = dolfinx.fem.Function(V)
f.interpolate(lambda x: (0.1*x[1]*abs(x[0]), 0.05*np.sin(2*x[0])))

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_vertices = dolfinx.mesh.compute_incident_entities(mesh.topology, boundary_facets, mesh.topology.dim-1, 0)
vertex_to_geometry = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object, 0, boundary_vertices, False)


c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
mesh.topology.create_connectivity(0, mesh.topology.dim)
v_to_c = mesh.topology.connectivity(0, mesh.topology.dim)

num_dofs = V.dofmap.index_map.size_local+ V.dofmap.index_map.num_ghosts
dof_to_geometry_map = np.full(num_dofs, -1, dtype=np.int32)
dofmap = V.dofmap
layout = dofmap.dof_layout

for i, (vertex, node) in enumerate(zip(boundary_vertices, vertex_to_geometry)):
    cell = v_to_c.links(vertex)[0]
    cell_dofs = dofmap.cell_dofs(cell)
    cvs = c_to_v.links(cell)
    local_index = np.flatnonzero(cvs == vertex)[0]
    dof = cell_dofs[layout.entity_dofs(0, local_index)]
    dof_to_geometry_map[dof] = node

# Create boundary perturbation vector
bs = V.dofmap.bs
perturbation_data = np.zeros((len(boundary_vertices), 3), dtype=np.float64)
geom = np.zeros(len(boundary_vertices), dtype=np.int32)
c = 0
for i, node in enumerate(dof_to_geometry_map):
    if node != -1:
        perturbation_data[c, :bs] = f.x.array[bs*i:bs*(i+1)]
        geom[c] = node
        c+=1
mesh.geometry.x[geom]+= perturbation_data
with dolfinx.io.XDMFFile(mesh.comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
1 Like