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

1 Like
16

I need the interpolation matrix across two different meshes in a monolithic fluid-structure interaction problem, which is exactly the same as you do here. Have you ever tried to do it with petsc4py.PETSc.Mat object in the parallel FEniCSx environment?

17

Do you mean creating directly a PETSc Matrix?

Anyway, since I need to solve hundreds of frequencies, I solve the domain in serial and parallelize the frequencies. I mean, every core is solving the entire domain for a single frequency. This makes the parallelism scaling basically perfect.

Then no, I am not using the parallel environment.

18

Thanks! It works with two scalar fields created by FunctionSpace. I am wondering how it cope with two vector fields by VectorFunctionSpace and further two FunctionSpace created by MixedElement?

In your method, it only returns an interpolation matrix I with dimension n_1\times n_2, where n_1 and n_2 are the numbers of grid points of two meshes. It dismatches the dimension for two function spaces which are not created by FunctionSpace(mesh, ("Lagrange", x)).

19

Yeah this is what I needed in the first place.
I will add support for vector fields and mixed elements.

20

I know I’m asking for a lot, but I would immensly appreciate it if you could update this function for dolfinx v0.8.0. I tried on my own but am just not able to do it. Currently, specifically, I cannot find any documentation on, nor can figure out on my own, what did dofmap.links() do before? msh.geometry.dofmap does not support it any longer, because dofmap is now just a np.ndarray…