# Extracting surface part of a Function

**URL:** <https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768>\
**Category:** General\
**Created:** [June 16, 2025, 1:10pm UTC](https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768 "2025-06-16T13:10:11Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Psirus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psirus/32/624_2.png) [@Psirus](https://fenicsproject.discourse.group/u/Psirus)\
**Post date:** [June 16, 2025, 1:10pm UTC](https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768/1 "2025-06-16T13:10:11Z")

</div>

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?

```python
import numpy as np
from mpi4py import MPI

import dolfinx

msh = dolfinx.mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]], [5, 5, 5])

fdim = msh.topology.dim - 1
msh.topology.create_connectivity(2, 3)
whole_exterior = dolfinx.mesh.exterior_facet_indices(msh.topology)

V = dolfinx.fem.functionspace(msh, ("CG", 1, (3,)))
u = dolfinx.fem.Function(V)

def stuff(x):
    return [np.sin(2.0*np.pi*x[0]), np.cos(2.0*np.pi*x[1]), x[2]]

u.interpolate(stuff)

xdmf_full = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "full.xdmf", "w")
xdmf_full.write_mesh(msh)
xdmf_full.write_function(u)
xdmf_full.close()

submesh, facet_map, vertex_map, _ = dolfinx.mesh.create_submesh(msh, fdim, whole_exterior)
V_sub = dolfinx.fem.functionspace(submesh, ("CG", 1, (3,)))
u_surf = dolfinx.fem.Function(V_sub)

parent_dofs = dolfinx.fem.locate_dofs_topological(V, fdim, whole_exterior)
indices = np.ravel([3*parent_dofs + i for i in range(3)], order="F")

u_surf.x.array[:] = u.x.array[indices]

xdmf_surf = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "surf.xdmf", "w")
xdmf_surf.write_mesh(submesh)
xdmf_surf.write_function(u_surf)

```

Y-component of full mesh on the left, my attempt of surface mesh on the right.

 ![Bildschirmfoto 2025-06-16 um 15.08.54](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/3/3e2468b243af5341f8ecc081f8fc2fef4afd2b57.jpeg)

---

<div class="post-metadata">

**Author:** ![Monisaddam](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/monisaddam/32/8536_2.png) [@Monisaddam](https://fenicsproject.discourse.group/u/Monisaddam)\
**Post date:** [June 16, 2025, 1:40pm UTC](https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768/2 "2025-06-16T13:40:29Z")

</div>

![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/574f0ea3c87b90159f7735c14762994e05de16bd.png)  
Hello, dear can you tell me how you are getting this cube??. I have save the file and open it in paraview but not getting your figure( Cube type).

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 16, 2025, 3:12pm UTC](https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768/3 "2025-06-16T15:12:32Z")

</div>

I’ve got a function for this, that hopefully will be merged into `scifem` eventually:

```python
"""
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()

```

---

<div class="post-metadata">

**Author:** ![Psirus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psirus/32/624_2.png) [@Psirus](https://fenicsproject.discourse.group/u/Psirus)\
**Post date:** [June 16, 2025, 3:29pm UTC](https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768/4 "2025-06-16T15:29:38Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 16, 2025, 4:03pm UTC](https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768/5 "2025-06-16T16:03:01Z")

</div>

> [@Psirus](#):
>
> `compute_entity_data`

Here is an implementation of that function

```python
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

```

---

<div class="post-metadata">

**Author:** ![Psirus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psirus/32/624_2.png) [@Psirus](https://fenicsproject.discourse.group/u/Psirus)\
**Post date:** [June 17, 2025, 7:33am UTC](https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768/6 "2025-06-17T07:33:29Z")

</div>

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.

```python
import numpy as np
from mpi4py import MPI

import dolfinx

from interpolate_submesh import interpolate_to_submesh

msh = dolfinx.mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]], [5, 5, 5])

fdim = msh.topology.dim - 1
msh.topology.create_connectivity(2, 3)
whole_exterior = dolfinx.mesh.exterior_facet_indices(msh.topology)

V = dolfinx.fem.functionspace(msh, ("CG", 1, (3,)))
u = dolfinx.fem.Function(V)

def stuff(x):
    return [np.sin(2.0*np.pi*x[0]), np.cos(2.0*np.pi*x[1]), x[2]]

u.interpolate(stuff)

xdmf_full = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "full.xdmf", "w")
xdmf_full.write_mesh(msh)
xdmf_full.write_function(u)
xdmf_full.close()

submesh, facet_map, vertex_map, _ = dolfinx.mesh.create_submesh(msh, fdim, whole_exterior)
V_sub = dolfinx.fem.functionspace(submesh, ("CG", 1, (3,)))
u_surf = dolfinx.fem.Function(V_sub)

interpolate_to_submesh(u, u_surf, facet_map)

xdmf_surf = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "surf.xdmf", "w")
xdmf_surf.write_mesh(submesh)
xdmf_surf.write_function(u_surf)

```

 ![Bildschirmfoto 2025-06-17 um 09.31.03](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/4/431100af5b0dcaeaee33a8d276630768e1bf326f.jpeg)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 22, 2025, 7:57pm UTC](https://fenicsproject.discourse.group/t/extracting-surface-part-of-a-function/17768/7 "2025-06-22T19:57:15Z")

</div>

This is due to a bug on v0.9.0, that has been fixed on the nightly branch.  
You can test this with for instance: `ghcr.io/fenics/dolfinx/dolfinx:nightly`