# Project submesh solution to original mesh

**URL:** <https://fenicsproject.discourse.group/t/project-submesh-solution-to-original-mesh/15340>\
**Category:** dolfinx\
**Created:** [July 25, 2024, 10:37am UTC](https://fenicsproject.discourse.group/t/project-submesh-solution-to-original-mesh/15340 "2024-07-25T10:37:29Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![hussleJ6](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/husslej6/32/4960_2.png) [@hussleJ6](https://fenicsproject.discourse.group/u/hussleJ6)\
**Post date:** [July 25, 2024, 10:37am UTC](https://fenicsproject.discourse.group/t/project-submesh-solution-to-original-mesh/15340/1 "2024-07-25T10:37:29Z")

</div>

Hi there,  
I had a question regarding the projection of a function/solution, from a submesh to the original mesh. It is comparable to a post before, so i will make use of this example;

> [@Extend submesh function to full mesh](https://fenicsproject.discourse.group/t/extend-submesh-function-to-full-mesh/8006):
>
> Hello ! I have a function defined on a submesh and I would like to extend it to the mesh that was used initially to create the submesh (with the value 0). Here’s a MWE of what I want to do (using main branch of dolfinx) import numpy as np from mpi4py import MPI import dolfinx, dolfinx.io, dolfinx.fem as fem, dolfinx.mesh mesh = dolfinx.mesh.create\_unit\_cube(MPI.COMM\_WORLD, 10, 10, 10, dolfinx.mesh.CellType.tetrahedron) submesh\_entities = dolfinx.mesh.locate\_entities(mesh, dim=3, marker=lamb…

In this discussion, a solution on a 3d submesh was returned into the original 3d mesh. Right now i want to accomplish to make a projection of a surface back into the original mesh.  
For the example provided, everything works fine;

```auto
import numpy as np
from mpi4py import MPI
import dolfinx
import dolfinx.io
import dolfinx.fem as fem
import dolfinx.mesh
import ufl

import pyvista
from dolfinx.plot import vtk_mesh

# Create a unit cube mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10, dolfinx.mesh.CellType.tetrahedron)

# Define a marker function to select the submesh domain
def marker(x):
    return (x[0] < 0.5) & (x[1] < 0.5)

# Locate entities based on the marker
submesh_entities = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, marker)

# Create submesh
submesh = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, submesh_entities)

# Function expression
def f_expr(x):
    return x

# If submesh returns a tuple, let's unpack it correctly:
if isinstance(submesh, tuple):
    submesh = submesh[0] # Assuming the first element is the mesh object

# Verify submesh type
print("Submesh type:", type(submesh))

element = ufl.VectorElement("CG", submesh.ufl_cell(), 1)
CG1_submesh = fem.FunctionSpace(submesh, element)

#CG1_submesh = fem.VectorFunctionSpace(submesh, ("CG", 1))
f_submesh = fem.Function(CG1_submesh)
f_submesh.interpolate(f_expr)

with dolfinx.io.XDMFFile(submesh.comm, "f_submesh2.xdmf", "w") as file:
    file.write_mesh(submesh)
    file.write_function(f_submesh)

element_mesh = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
CG1_mesh = fem.FunctionSpace(mesh, element_mesh)

#CG1_mesh = fem.VectorFunctionSpace(mesh, ("CG", 1))
f_mesh = fem.Function(CG1_mesh)
f_mesh.vector.set(0.)

# Setup for interpolation
tree = dolfinx.geometry.bb_tree(submesh, 3)
num_entities_local = submesh.topology.index_map(3).size_local + submesh.topology.index_map(3).num_ghosts
entities = np.arange(num_entities_local, dtype=np.int32)
midpoint_tree = dolfinx.geometry.create_midpoint_tree(submesh, 3, entities)

def f_mesh_expr(x):
    cells = dolfinx.geometry.compute_closest_entity(tree, midpoint_tree, submesh, x.T)
    return np.ascontiguousarray(f_submesh.eval(x.T, cells).T)

f_mesh.interpolate(f_mesh_expr, cells=submesh_entities)

with dolfinx.io.XDMFFile(mesh.comm, "f_mesh.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(f_mesh)

```

If I try to do this for a surface, it seems to write the surface correctly to a .xdmf file. Still it is not able to project the solution back to the original mesh. Likely this has to do with some surface type of input where it would expect a cell, but I’m not able to figure out how to adjust this. Any suggestion would be appreciated.

```auto
import numpy as np
from mpi4py import MPI
import dolfinx
import dolfinx.io
import dolfinx.fem as fem
import dolfinx.mesh
import ufl

# Create a unit cube mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10, dolfinx.mesh.CellType.tetrahedron)

# Define a marker function to select the submesh domain
def boundary_marker(x):
    # x[1] represents the y-coordinate of the midpoint of each facet
    return x[1] < 0.1

# Locate boundary entities
fdim = mesh.topology.dim - 1 # Dimension of the facets
boundary_entities = dolfinx.mesh.locate_entities_boundary(mesh, fdim, boundary_marker)
# Ensure connectivity is built for the facets to be extracted
mesh.topology.create_connectivity(fdim, mesh.topology.dim)

# Create submesh
boundary_mesh = dolfinx.mesh.create_submesh(mesh, fdim, boundary_entities)

# Function expression
def f_expr(x):
    return x

# If submesh returns a tuple, let's unpack it correctly:
if isinstance(boundary_mesh, tuple):
    boundary_mesh = boundary_mesh[0] # Assuming the first element is the mesh object
    
# Check the type of the boundary mesh and print its topology dimension
print("Boundary mesh topology dimension:", boundary_mesh.topology.dim)

# Verify submesh type
print("Submesh type:", type(boundary_mesh))

element = ufl.VectorElement("CG", boundary_mesh.ufl_cell(), 1)
CG1_submesh = fem.FunctionSpace(boundary_mesh, element)

#CG1_submesh = fem.VectorFunctionSpace(submesh, ("CG", 1))
f_submesh = fem.Function(CG1_submesh)
f_submesh.interpolate(f_expr)

# Writing the boundary mesh to a file
with dolfinx.io.XDMFFile(boundary_mesh.comm, "boundary_mesh.xdmf", "w") as file:
    file.write_mesh(boundary_mesh)
    file.write_function(f_submesh)

element_mesh = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
CG1_mesh = fem.FunctionSpace(mesh, element_mesh)

#CG1_mesh = fem.VectorFunctionSpace(mesh, ("CG", 1))
f_mesh = fem.Function(CG1_mesh)
f_mesh.vector.set(0.)

tree = dolfinx.geometry.bb_tree(boundary_mesh, 2)
num_entities_local = boundary_mesh.topology.index_map(2).size_local + boundary_mesh.topology.index_map(2).num_ghosts
entities = np.arange(num_entities_local, dtype=np.int32)
midpoint_tree = dolfinx.geometry.create_midpoint_tree(boundary_mesh, 2, entities)

def f_mesh_expr(x):
    cells = dolfinx.geometry.compute_closest_entity(tree, midpoint_tree, boundary_mesh, x.T)
    return np.ascontiguousarray(f_submesh.eval(x.T, cells).T)

f_mesh.interpolate(f_mesh_expr, cells=boundary_entities)

with dolfinx.io.XDMFFile(mesh.comm, "f_mesh.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(f_mesh)

```

---

<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:** [August 23, 2024, 6:40am UTC](https://fenicsproject.discourse.group/t/project-submesh-solution-to-original-mesh/15340/2 "2024-08-23T06:40:30Z")

</div>

I do not see any projection operator in this code. Please elaborate on what you classify as a projection operator. You should also consider

> <https://github.com/FEniCS/dolfinx/pull/3114>
>
> A feature that is quite neat for post-processing, as one often only wants to vis…ualize on a sub mesh.
> Can also be used to avoid mixed-assembly in cases of iterative schemes:
> 1. Solve on a parent mesh
> 2. Move solution to sub mesh
> 3. Solve equation on sub mesh
> 4. Move solution back and repeat.
> 
> There is a slight API change in the C++ interface, as expressions need the associated mesh to be supplied.
> In Python, the API only changes if you want to map between meshes.

and the subsequent refactoring

> <https://github.com/FEniCS/dolfinx/pull/3177>
>
> Resolves #3173.
> Simplifies logic of interpolate.
> \*\*API Change\*\*
> Non-matching …meshes interpolation now has its own interface: \`interpolate\_nonmatching\_meshes\` in both C++ and Python.
> It also clarifies that interpolation of \`dolfinx::fem::Expression\` is not supported on non-matching meshes.