I have a vector-valued function u defined on the full mesh. Now I want to export just (a part of) the surface of the mesh, including this function u. My attempt so far looks like this, but the result isn’t correct. Do I need some sort of reordering in the LHS of u_sarf.x.array[:] = ... instead of :? Or is indices = np.ravel([3*parent_dofs + i for i in range(3)], order="F")not correct? BTW, what is the recommended way to go from “dofs” to actual array indices for vector-valued functions?
I’ve got a function for this, that hopefully will be merged into scifem eventually:
"""
Interpolate from a parent mesh to a set of exterior facets
Author: Jørgen S. Dokken
SPDX License identifier: MIT
"""
def interpolate_to_submesh(
u_volume: dolfinx.fem.Function,
u_surface: dolfinx.fem.Function,
entity_map: npt.NDArray[np.int32],
):
"""
Interpolate a function `u_volume` into the function `u_surface`.
Note:
Does not work for DG as no dofs are associated with the facets
Args:
u_volume: Function to interpolate data from
u_surface: Function to interpolate data to
entity_map: Map from cells in the surface mesh to facets in the volume mesh
"""
V_vol = u_volume.function_space
mesh = V_vol.mesh
V_surf = u_surface.function_space
submesh = V_surf.mesh
if Version(dolfinx.__version__) > Version("0.9.0"):
ip = V_surf.element.interpolation_points
else:
ip = V_surf.element.interpolation_points()
expr = dolfinx.fem.Expression(u_volume, ip)
mesh.topology.create_connectivity(mesh.topology.dim, submesh.topology.dim)
mesh.topology.create_connectivity(submesh.topology.dim, mesh.topology.dim)
c_to_e = mesh.topology.connectivity(mesh.topology.dim, submesh.topology.dim)
e_to_c = mesh.topology.connectivity(submesh.topology.dim, mesh.topology.dim)
# Pack integration entities, including ghosts
integration_entities = np.zeros((len(entity_map), 2), dtype=np.int32)
for i, facet in enumerate(entity_map):
cells = e_to_c.links(facet)
if len(cells) > 1:
raise RuntimeError("Only exterior facets are supported for now")
facets = c_to_e.links(cells[0])
local_index = np.flatnonzero(facets == facet)
assert len(local_index) == 1
integration_entities[i] = (cells[0], local_index[0])
# Data is (num_cells, num_points/num_dofs, value_size)
if Version(dolfinx.__version__) > Version("0.9.0"):
data = expr.eval(mesh, integration_entities)
else:
data = expr.eval(mesh, integration_entities.flatten())
submesh.topology.create_entity_permutations()
mesh.topology.create_entity_permutations()
cell_info = mesh.topology.get_cell_permutation_info()
ft = V_surf.element.basix_element.cell_type
for i in range(integration_entities.shape[0]):
perm = np.arange(data.shape[1], dtype=np.int32)
try:
# Can be used once https://github.com/FEniCS/basix/pull/913
V_vol.element.basix_element.permute_subentity_closure_inv(
perm, cell_info[integration_entities[i, 0]], ft, int(integration_entities[i, 1])
)
except TypeError:
entity_info = compute_entity_data(
V_vol.element, cell_info[integration_entities[i, 0]], ft, integration_entities[i, 1]
)
V_vol.element.basix_element.permute_subentity_closure_inv(perm, entity_info, ft)
data[i] = data[i][perm]
if len(data.shape) == 3:
# Data is now (num_cells, value_size,num_points)
data = data.swapaxes(1, 2)
# Data is now (value_size, num_cells, num_points)
data = data.swapaxes(0, 1)
if expr.value_size == 1:
shaped_data = data.flatten()
else:
shaped_data = data.reshape(expr.value_size, -1)
u_surface._cpp_object.interpolate(shaped_data, np.arange(len(entity_map), dtype=np.int32))
u_surface.x.scatter_forward()
Thank you for taking the time, this is much more involved than I realised. For me it triggers the TypeError, which is expected I guess, but I don’t know where compute_entity_data is from, so I can’t import it.
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem
import dolfinx
import numpy.typing as npt
import basix
from packaging.version import Version
def compute_entity_data(
element,
cell_info: int,
entity_type: basix.CellType,
entity_index: int,
) -> int:
"""
Compute entity permutation data on a sub entity of the mesh.
Note:
This can be removed once https://github.com/FEniCS/basix/pull/913 is merged.
Args:
element: The element to compute permutations for
cell_info: Permutation data for the cell
entity_type: The type of the entity
entity_index: The local index of the entity (relative to the cell)
Returns:
entity_info: Permutation data for the entity
"""
ct = element.basix_element.cell_type
tdim = len(basix.cell.topology(ct)) - 1
num_entity_dofs = len(element.basix_element.entity_dofs[2])
face_start = 3 * num_entity_dofs if tdim == 3 else 0
entity_info = 0
edim = len(basix.cell.topology(entity_type)) - 1
if edim == 0:
pass
elif edim == 1:
entity_info = cell_info >> (face_start + entity_index) & 1
elif edim == 2:
entity_info = cell_info >> (3 * entity_index) & 7
else:
raise ValueError(f"Unsupported entity {entity_type}")
return entity_info
It runs now, but still doesn’t give me the correct result. I’m assuming the entity map I have to pass in is the entity map returned by create_submesh, which I’ve called facet_map here. Is that correct? I’m using dolfinx 0.9.0.
