Map function on boundary submesh to full mesh

dolfinx
1

Hello,

I have the following problem. I have some finite element problem where the weak formulation contains an integral over the domain boundary. This integral contains a function that is defined on the boundary. Say something that looks like,

\int_{\partial \Omega} dS \, (\vec{n} \cdot \vec{\nabla} u) v f(\vec{x})

Where u and v are the trial and test functions respectively and f: \partial \Omega \to \mathbb{R}. The question is then how should one implement this. The following snippet of code does the job:

from mpi4py import MPI
import dolfinx.mesh as mesh
import dolfinx.fem as fem
import numpy as np
import numpy.typing as npt
import ufl

if __name__ == "__main__":
    # Define mesh
    n_cells: int = 20
    domain: mesh.Mesh = mesh.create_rectangle(
        comm=MPI.COMM_WORLD,
        points=[[-1.0, -1.0], [1.0, 1.0]],
        n=[n_cells, n_cells],
        dtype=np.float64,
    )
    tdim: int = domain.topology.dim

    # Create boundary mesh
    domain.topology.create_connectivity(tdim - 1, tdim)
    bndry_facets: npt.NDArray[np.int32] = mesh.exterior_facet_indices(domain.topology)
    bndry_domain, entity_map, vertex_map, node_map = mesh.create_submesh(
        msh=domain, dim=tdim - 1, entities=bndry_facets
    )
    bndry_facet_tags: mesh.MeshTags = mesh.meshtags(
        msh=domain,
        dim=tdim - 1,
        entities=bndry_facets,
        values=np.ones(bndry_facets.shape[0], dtype=np.int32),
    )

    # Define function spaces
    element: fem.ElementMetaData = fem.ElementMetaData(family="Lagrange", degree=1)
    V_full: fem.FunctionSpace = fem.functionspace(mesh=domain, element=element)
    V_bndry: fem.FunctionSpace = fem.functionspace(mesh=bndry_domain, element=element)

    fnc_bndry = fem.Function(V_bndry)
    fnc_full = fem.Function(V_full)
    assert isinstance(fnc_bndry, fem.Function) and isinstance(fnc_full, fem.Function)

    # Define function on the boundary
    def fnc_eval(x: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
        return np.sum(x, axis=0)

    fnc_bndry.interpolate(u0=fnc_eval)

    ### THIS IS THE IMPORTANT PART ###
    # Map function on the boundary to the entire mesh
    bndry_dofs: npt.NDArray[np.int32] = fem.locate_dofs_topological(
        V=V_full, entity_dim=tdim - 1, entities=bndry_facets
    )
    cell_map = domain.topology.index_map(tdim)
    num_cells: int = cell_map.size_local + cell_map.num_ghosts
    cells: npt.NDArray[np.int32] = np.arange(num_cells, dtype=np.int32)
    interpolation_data = fem.create_interpolation_data(
        V_to=V_full, V_from=V_bndry, cells=cells
    )
    fnc_full.interpolate_nonmatching(
        u0=fnc_bndry, cells=cells, interpolation_data=interpolation_data
    )
    ##################################

    # Weak formulation with boundary integral
    u = ufl.TrialFunction(V_full)
    v = ufl.TestFunction(V_full)
    n = ufl.FacetNormal(domain=domain)
    ds = ufl.Measure("ds", domain=domain, subdomain_data=bndry_facet_tags)

    some_ufl = fnc_full * ufl.inner(ufl.inner(n, ufl.grad(u)), v) * ds(subdomain_id=1)

We create sub-mesh, define f on this sub-mesh and then use the interpolate_nonmatching method to “lift” the function f from \partial \Omega to \Omega. This approach does not make use of the entity_map, vetex_map or node_map that are provided when creating the sub-mesh. I therefore wonder if there is a more efficient approach than then one presented here.

Note: The integration shown here is over the entire boundary. That may not always be the case. I would like an alternative solution to still work in those cases.

2

The problem seems very similar to the one discussed in this thread . If true, you do not need to lift the function from sub to main mesh in order to add the term in the weak form. Passing the entity maps to the form compiler should be enough as explained in the referred thread.

3

Thanks for the response, the answer was shown clearly in:

https://jsdokken.com/FEniCS-workshop/src/multiphysics/coupling.html