Extend submesh function to full mesh

dolfinx
1

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=lambda x: (x[0] < 0.5) & (x[1] < 0.5))
submesh, vertex_map, geom_map = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, submesh_entities)

def f_expr(x):
    return x[0]

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

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

tree = dolfinx.geometry.BoundingBoxTree(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 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)

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

This code actually works and we can see the results in paraview (submesh solution on the left and extended solution on the right).

mwe_1
mwe_11094×758 89.3 KB

However it fails if I try to do the same with vector valued functions :

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=lambda x: (x[0] < 0.5) & (x[1] < 0.5))
submesh, vertex_map, geom_map = dolfinx.mesh.create_submesh(mesh, mesh.topology.dim, submesh_entities)

def f_expr(x):
    return x

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

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

tree = dolfinx.geometry.BoundingBoxTree(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 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)

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

We can see that the values are not correctly computed

mwe_2
mwe_21031×772 128 KB

So obviously I must be doing something wrong but I can’t figure out what that would be.
Any help would be appreciated!
Thanks

2 Likes
2

This is actually really odd. I’m not yet certain why this is happening.

If you supply:

def f_mesh_expr(x):
    cells = dolfinx.geometry.compute_closest_entity(tree, midpoint_tree, submesh, x.T)
    u = f_submesh.eval(x.T, cells).T
    print("u == f_expr(x) :", np.all(np.isclose(f_expr(x), u)))
    return u

and interchange returning the equivalent arrays u and f_expr(x) it either works or doesn’t… bizarre…

3

I can file a bug on the github repo if necessary?
As a workaround I ended up doing this

def f_mesh_expr(x):
    cells = dolfinx.geometry.compute_closest_entity(tree, midpoint_tree, submesh, x.T)

    val = np.zeros_like(x)
    for i in range(np.shape(x)[1]):
        val[:, i] = f_submesh.eval(x[:, i], cells[i])
    return val

which works fine but has suboptimal performance since it uses a python loop.

1 Like
4

thanks for this. Let me see if I can figure out what’s going on with your update.

5

I think I’ve found the problem:

    print(val.flags)
    print(u.flags)
    print(x.flags)

gives us

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

Once the data is being passed back, if it’s not correctly aligned in terms of contiguity, then when it gets flattened in C++ the ordering will be incorrect.

I’ll raise an issue after figuring out at which point the C/F_CONTIGUOUS issue occurs.

Nice find.

4 Likes
6

So as a quick fix for your problem:

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)
1 Like
7

That was indeed quite a subtle bug, thanks for the fix!