Formulation with functions from parent mesh and submesh

dolfinx
1

Hi all,

I have a function defined on a submesh and a function defined on the parent mesh.

I would like to use both in a formulation (consider the following MWE).

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import ufl
import numpy as np
import basix


def half(x):
    return x[1] <= 0.5 + 1e-14


mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 10, 10, dolfinx.mesh.CellType.triangle
)

# Split domain in half
gdim = mesh.geometry.dim
tdim = mesh.topology.dim

bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half)
num_cells_local = (
    mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
)
cells = np.full(num_cells_local, 4, dtype=np.int32)
cells[bottom_cells] = 3

ct = dolfinx.mesh.meshtags(
    mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells
)

submesh_b, submesh_b_to_mesh, b_v_map = dolfinx.mesh.create_submesh(
    mesh, tdim, ct.find(3)
)[0:3]
submesh_t, submesh_t_to_mesh, t_v_map = dolfinx.mesh.create_submesh(
    mesh, tdim, ct.find(4)
)[0:3]

def create_function(mesh):
    degree = 1
    element_CG = basix.ufl.element(
        basix.ElementFamily.P,
        mesh.basix_cell(),
        degree,
        basix.LagrangeVariant.equispaced,
    )
    V = dolfinx.fem.functionspace(mesh, element_CG)
    u = dolfinx.fem.Function(V)
    return u

u_b = create_function(submesh_b)
u_t = create_function(submesh_t)

u_b.interpolate(lambda x: x[1])
u_t.interpolate(lambda x: 1 - x[1])

# Create a function on the full mesh
T = create_function(mesh)
T.interpolate(lambda x: 10 + x[0])

# compute something with a mix of both
ds_b = ufl.Measure("ds", domain=submesh_b)
dx_b = ufl.Measure("dx", domain=submesh_b)
dx = ufl.Measure("dx", domain=mesh)

n_b = ufl.FacetNormal(submesh_b)

form = dolfinx.fem.form(T*ufl.dot(ufl.grad(u_b), n_b) * ds_b)
print(dolfinx.fem.assemble_scalar(form))

This produces:

Traceback (most recent call last):
  File "/workspaces/FESTIM/mwe_post_processing.py", line 89, in <module>
    print(dolfinx.fem.assemble_scalar(form))
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/assemble.py", line 140, in assemble_scalar
    coeffs = coeffs or _pack_coefficients(M._cpp_object)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IndexError: map::at

I found that creating another function defined on the submesh and interpolating works:

T_b = dolfinx.fem.Function(u_b.function_space)
T_b.interpolate(T)
form = dolfinx.fem.form(T_b*ufl.dot(ufl.grad(u_b), n_b) * ds_b)
print(dolfinx.fem.assemble_scalar(form))

I was wondering if this was the way to go or if there was another (better?) method.

Cheers!

2

You need to supply entity maps from the integration domain, here

The domains submesh_b. Therefore you need to send in the sub_to_parent_entity map, ie.

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import ufl
import numpy as np
import basix


def half(x):
    return x[1] <= 0.5 + 1e-14


mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 10, 10, dolfinx.mesh.CellType.triangle
)

# Split domain in half
gdim = mesh.geometry.dim
tdim = mesh.topology.dim

bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half)
num_cells_local = (
    mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
)
cells = np.full(num_cells_local, 4, dtype=np.int32)
cells[bottom_cells] = 3

ct = dolfinx.mesh.meshtags(
    mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells
)

submesh_b, submesh_b_to_mesh, b_v_map = dolfinx.mesh.create_submesh(
    mesh, tdim, ct.find(3)
)[0:3]


def create_function(mesh):
    degree = 1
    element_CG = basix.ufl.element(
        basix.ElementFamily.P,
        mesh.basix_cell(),
        degree,
        basix.LagrangeVariant.equispaced,
    )
    V = dolfinx.fem.functionspace(mesh, element_CG)
    u = dolfinx.fem.Function(V)
    return u


u_b = create_function(submesh_b)

u_b.interpolate(lambda x: x[1])

# Create a function on the full mesh
T = create_function(mesh)
T.interpolate(lambda x: 10 + x[0])
T.name = "T"
# compute something with a mix of both
ds_b = ufl.Measure("ds", domain=submesh_b)
dx_b = ufl.Measure("dx", domain=submesh_b)
dx = ufl.Measure("dx", domain=mesh)

n_b = ufl.FacetNormal(submesh_b)
u_b.name = "u_b"
form = dolfinx.fem.form(
    T * ufl.dot(ufl.grad(u_b), n_b) * ds_b,
    entity_maps={mesh: submesh_b_to_mesh},
)
print(dolfinx.fem.assemble_scalar(form))
3

Which version of dolfinx is required to execute this code? I just tried on 0.8.0 (conda) and jpdean’s branch on docker (which is at 0.7.0) but none worked.

4

The dolfinx nightly docker image works, or the 0.8.0 image, if you change:

to

form = dolfinx.fem.form(
    T * ufl.dot(ufl.grad(u_b), n_b) * ds_b,
    entity_maps={mesh._cpp_object: submesh_b_to_mesh},
)
1 Like