Project submesh solution to original mesh

dolfinx
1

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;

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;

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.

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)
2

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

and the subsequent refactoring