Interpolation matrix with non matching meshes

Linear Algebra
1

Hello guys!

Since it is possible to interpolate one function defined on a function space into another one defined on another function space (from a non matching mesh), like so:

V_1 = dolfinx.fem.FunctionSpace(msh_1, ("CG", 1))
V_2 = dolfinx.fem.FunctionSpace(msh_2, ("CG", 1))

u_1 = dolfinx.fem.Function(V_1)
u_2 = dolfinx.fem.Function(V_2)

#... compute or assign u_1 ...

u_2.interpolate(u_1)

I wonder if I can use the interpolate function to build an interpolation matrix that links the two function vectors u_1.vector and u_2.vector such that

\bf{u_2} = \bf{I} \bf{u_1}

where \bf{I} Is a rectangular matrix, such that \bf{u_2} is the same that I would obtain from the interpolation function

Thank you very much

2

See https://github.com/FEniCS/dolfinx/blob/e4439ccca81b976d11c6f606d9c612afcf010a31/python/test/unit/fem/test_petsc_discrete_operators.py#L134

3

Looks cool! I’ll try it!

4

I am getting

G = interpolation_matrix(V_a._cpp_object, V_s._cpp_object)
RuntimeError: Cannot insert rows that do not exist in the IndexMap.
5

Oh, sorry, i didn’t see that your are doing interpolation across different meshes.

There is currently not a way or creating that matrix without doing it by hand, but you can cache alot of the data/communication required for the cross mesh interpolation, see:

6

Ok, thank you Jørgen, I’ll have a look.

P.S. Are these new features of 0.7.0 version?

7

I think the non-matching interpolation was introduced in 0.6.0 (@Massimiliano_Leoni can comment on that). The caching was introduced in 0.7 I think.

8

create_nonmatching_meshes_interpolation_data() seems to associate the point where the function has to be interpolated with the cell that contains it. It is indeed not varying when the two meshes don’t change and saves a lot of computations for multiple interpolations.

Is there a way that allows me, given a cell and a point contained in it, to retrieve the shape functions and the corresponding cell nodes?

9

You can check out how dolfinx.fem.Function.eval works, as it does all the steps needed to pull back points from physical space to the reference, and use the underlying dolfinx.fem.FiniteElement to tabulate the basis functions as these points.

10

Thank you again,
I will try it!!

11

Hello Jørgen,

I am almost done, but I am struggling with the function msh.geometry.cmap.pull_back

I see that you wrote some working code:

That works fine with one point. With multiple points, the function seems not working. I think that its inputs shape matters.

In particular, how can I modify the following:

cell = cells.links(0)[0]

cmap = domain.geometry.cmap
geom_dofs = domain.geometry.dofmap.links(cell)
x_ref = cmap.pull_back(point, domain.geometry.x[geom_dofs])

when cell is an array of multiple cells?

Thank you very much,

Antonio

12

It only supports pull back on single cells (but for multiple points in that cell): https://github.com/FEniCS/dolfinx/blob/26fbb085d5d3c39340dab91705e4c475e98f4ec5/python/dolfinx/wrappers/fem.cpp#L737-L759

13

Yeah, now the errors make sense, thank you.

The working slice of code for multiple points is the following:

for i in range(0, len(cells)):
    geom_dofs = msh_a.geometry.dofmap.links(cells[i])
    x_ref[i,:] = msh_a.geometry.cmap.pull_back([points[i,:]], x[geom_dofs])

Sorry for the lack of context, but after I clean the complete code for the interpolation matrix, I am going to put it here.

14

The interpolation matrix function for nonmatching meshes is done!
Works on dolfinx v0.6.0
Works only for lagrange element

import numpy as np
import dolfinx
from dolfinx import geometry
import basix


def interpolation_matrix_nonmatching_meshes(V_1,V_0): # Function spaces from nonmatching meshes
    msh_0 = V_0.mesh
    msh_1 = V_1.mesh
    x_0   = V_0.tabulate_dof_coordinates()
    x_1   = V_1.tabulate_dof_coordinates()

    bb_tree         = geometry.BoundingBoxTree(msh_0, msh_0.topology.dim)
    cell_candidates = geometry.compute_collisions(bb_tree, x_1)
    cells           = []
    points_on_proc  = []
    index_points    = []
    colliding_cells = geometry.compute_colliding_cells(msh_0, cell_candidates, x_1)

    for i, point in enumerate(x_1):
        if len(colliding_cells.links(i))>0:
            points_on_proc.append(point)
            cells.append(colliding_cells.links(i)[0])
            index_points.append(i)
            
    index_points_   = np.array(index_points)
    points_on_proc_ = np.array(points_on_proc, dtype=np.float64)
    cells_          = np.array(cells)

    ct      = dolfinx.cpp.mesh.to_string(msh_0.topology.cell_type)
    element = basix.create_element(basix.finite_element.string_to_family(
        "Lagrange", ct), basix.cell.string_to_type(ct), V_0.ufl_element().degree(), basix.LagrangeVariant.equispaced)

    x_ref = np.zeros((len(cells_), 3))

    for i in range(0, len(cells_)):
        geom_dofs  = msh_0.geometry.dofmap.links(cells_[i])
        x_ref[i,:] = msh_0.geometry.cmap.pull_back([points_on_proc_[i,:]], msh_0.geometry.x[geom_dofs])

    basis_matrix = element.tabulate(0, x_ref)[0,:,:,0]

    cell_dofs         = np.zeros((len(x_1), len(basis_matrix[0,:])))
    basis_matrix_full = np.zeros((len(x_1), len(basis_matrix[0,:])))


    for nn in range(0,len(cells_)):
        cell_dofs[index_points_[nn],:] = V_0.dofmap.cell_dofs(cells_[nn])
        basis_matrix_full[index_points_[nn],:] = basis_matrix[nn,:]

    cell_dofs_ = cell_dofs.astype(int) ###### REDUCE HERE

    I = np.zeros((len(x_1), len(x_0)), dtype=complex)

    for i in range(0,len(x_1)):
        for j in range(0,len(basis_matrix[0,:])):
            I[i,cell_dofs_[i,j]] = basis_matrix_full[i,j]

    return I

You can find it also here:

15

My function works the same of the following:

but here the FunctionSpaces come from two different meshes.
Suppose that we want to interpolate a function from V_0 to V_1. It results:

\mathbf{I} = f(V_1, V_0)

Now, suppose that we have the following situation:

V_0 = FunctionSpace(mesh_0,("CG",deg_0))
V_1 = FunctionSpace(mesh_1,("CG",deg_1))

p_0 = Function(V_0)
...
p_0 = SomeComputations()

we define a new function on V_1:

p_1 = Function(V_1)

If we wanted to interpolate the function, we could use the command:

p_1.interpolate(p_0)

But now we have an alternative. We can indeed compute the interpolation matrix and compute the p_1 array from a matrix multiplication between I_matrix and the p_0 array:

I_matrix = interpolation_matrix_nonmatching_meshes(V_1,V_0)
p_1.x.array[:] = np.matmul(I_matrix,p_0.x.array)

getting the same result.

The purpose of such matrix is to try building a monolithic coupled vibro-acoustic system where the domains meshes at the interface don’t match.

Maybe someone else finds this useful.

Bye guys, I need a beer now.

Antonio