Interpolate function from lower to higher dimension mesh

dolfinx
1

Hi,

I need to interpolate a function defined on a submesh to the main mesh of higher dimension. I can make it work for nodal elements but I need to get it to work for functions defined on Nedelec space. I tried to use the functions introduced in release 0.9 but there seems to be an error I do not understand. I would appreciate if someone could fix the following MWE:

import dolfinx, ufl, gmsh, mpi4py
import numpy as np
from dolfinx.io import gmshio

def create_geometry():
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0) # disable output messages
    gmsh.clear()
    gdim = 2
    model_rank = 0

    occ = gmsh.model.occ
    gdim = 2

    fov = occ.add_rectangle(0, 0, 0, 1, 1)
    
    occ.synchronize()
    # number all domains
    all_doms = gmsh.model.getEntities(gdim)
    for j, dom in enumerate(all_doms):
        gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

    # number all boundaries
    all_edges = gmsh.model.getEntities(gdim - 1)
    for j, edge in enumerate(all_edges):
        gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

    # gmsh.model.mesh.setSize(gmsh.model.getEntities(), max_size)
    gmsh.option.setNumber("Mesh.MeshSizeMax", 0.1)

    gmsh.model.mesh.generate(gdim)
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
    return mesh, ct, ft

mesh, ct, ft = create_geometry()

# create submesh at x = 0
marker = lambda x: np.isclose(x[0], 0.0)
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, marker)
submesh, entity_map, vertex_map, node_map = dolfinx.mesh.create_submesh(mesh, 1, left_facets)

V_submesh = dolfinx.fem.functionspace(submesh, ("CG", 2, (2,))) # a vector function
uh_submesh = dolfinx.fem.Function(V_submesh)

V = dolfinx.fem.functionspace(mesh, ("N1curl", 2))
uh = dolfinx.fem.Function(V)

submesh_to_mesh_interp_data = dolfinx.fem.create_interpolation_data(V, V_submesh, ct.find(1)) # I only care about the boundary so take all cells for simplicity

# TOFIX:
uh.interpolate_nonmatching(uh_submesh, cells=None, interpolation_data=submesh_to_mesh_interp_data)
3

Why do you need to interpolate the solution back to the parent grid?

You can use function from the submesh in a form with the parent mesh.

4

Could you please point out the syntax of adding this boundary field to the weak form of the main mesh? This is indeed what I want to do for now.

5

See for instance
http://jsdokken.com/FEniCS-workshop/src/multiphysics/coupling.html

6

Thanks for the link. I adapted it into my code but it leads to error when I solve the problem. The following MWE reproduces the problem (conda, 0.9, complex):

import dolfinx, ufl, gmsh, mpi4py
import numpy as np
from dolfinx.io import gmshio
from dolfinx.fem.petsc import LinearProblem


def create_geometry():
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0) # disable output messages
    gmsh.clear()
    gdim = 2
    model_rank = 0

    occ = gmsh.model.occ
    gdim = 2

    fov = occ.add_rectangle(0, 0, 0, 1, 1)
    
    occ.synchronize()
    # number all domains
    all_doms = gmsh.model.getEntities(gdim)
    for j, dom in enumerate(all_doms):
        gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

    # number all boundaries
    all_edges = gmsh.model.getEntities(gdim - 1)
    for j, edge in enumerate(all_edges):
        gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

    # gmsh.model.mesh.setSize(gmsh.model.getEntities(), max_size)
    gmsh.option.setNumber("Mesh.MeshSizeMax", 0.1)

    gmsh.model.mesh.generate(gdim)
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
    return mesh, ct, ft

mesh, ct, ft = create_geometry()
submesh_bnd_idx = 4 # boundary where submesh is created at x=0, found by inspection

# create submesh at x = 0
marker = lambda x: np.isclose(x[0], 0.0)
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, marker)
submesh, entity_map, vertex_map, node_map = dolfinx.mesh.create_submesh(mesh, 1, left_facets)

V_submesh = dolfinx.fem.functionspace(submesh, ("CG", 2)) # a vector function
uh_submesh = dolfinx.fem.Function(V_submesh)

V = dolfinx.fem.functionspace(mesh, ("CG", 2))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# submesh_to_mesh_interp_data = dolfinx.fem.create_interpolation_data(V, V_submesh, ct.find(1))
# # TOFIX:
# uh.interpolate_nonmatching(uh_submesh, cells=None, interpolation_data=submesh_to_mesh_interp_data)

facet_imap = mesh.topology.index_map(ft.dim)
num_facets = facet_imap.size_local + facet_imap.num_ghosts
mesh_to_submesh = np.full(num_facets, -1)
mesh_to_submesh[entity_map] = np.arange(len(entity_map))
entity_maps = {submesh: mesh_to_submesh}

ds_bnd = ufl.Measure("ds", domain=mesh, subdomain_data=ft, subdomain_id=submesh_bnd_idx)

qs = dolfinx.fem.Constant(mesh, 0.0j)
F = -ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx + ufl.inner(u, v)*ufl.dx + ufl.inner(qs, v)*ufl.dx
F += ufl.inner(uh_submesh, v) * ds_bnd

a = ufl.lhs(F)
L = ufl.rhs(F)

a_form = dolfinx.fem.form(a)
L_form = dolfinx.fem.form(L, entity_maps=entity_maps)

problem = LinearProblem(a_form, L_form, bcs=[])
uh = problem.solve()
7

I realized that when passing form object to LinearProblem, the solution function needs to be specified as well. It works now, thanks! The following correction to the code above fixes it:

uh = dolfinx.fem.Function(V)
problem = LinearProblem(a_form, L_form, u=uh, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})

It is unclear to me though how this approach will extend if I have two submeshes? How will I pass entity_maps to the form compiler?

8

You pass one entity map per sub-mesh. i.e.
entity_maps = {submesh0: mesh_to_submesh0, submesh1: mesh_to_submesh1}

1 Like