High order interpolation between two meshes

Hello,
I’m working on a parallel-in-time method and for the spatial discretization I’m using dolfinx. The method that I use is multilevel and employ a hierarchy of meshes in a multigrid fashion (PFASST). It is known that for efficiency the prolongation interpolator should be high order. So I was wondering how can be interpolate a function lying on a coarse mesh onto a fine mesh but using high order interpolation?

Currently I’m doing something like

u_fine.interpolate(u_coarse,
                   nmm_interpolation_data=fem.create_nonmatching_meshes_interpolation_data(
                                     u_fine.function_space.mesh._cpp_object, 
                                     u_fine.function_space.element, 
                                     u_coarse.function_space.mesh._cpp_object))

In my case both functions u_fine and u_coarse are on a space with Lagrangian basis functions of order 1. I wanted to ask which is the order of the inteprolation that I’m doing? One? Is it possible to increase it?

Thank you for your time :slight_smile:

What do you imply by the efficiency of the interpolator?

Do mean computational speed or accuracy?

Sorry, I meant efficiency of the parallel in time method. This method is similar to multigrid and does V-cycles. However, when I go from a coarse mesh to fine one the inteprolation must be high order or the number of iterations increase significantly.

I still don’t quite get your question.

Interpolation itself scales with the number of degrees of freedom (linearly).
The non-matching interpolation is split into two steps:

  1. For all interpolation points on your find grid find what processes on the coarse grid owns these points. This is what fem.create_nonmatching_meshes_interpolation_data( u_fine.function_space.mesh._cpp_object, u_fine.function_space.element, u_coarse.function_space.mesh._cpp_object) does. If your meshes are not moving in time, then this function can be called once.
  2. Each process evaluates the coarse function at the points that is assign to it. This process involves pulling back the physical points to the reference cell (which involves a Newton solve for non-affine geometries) and tabulating and combining basis functions. Then these values are sent to the relevant processes using MPI.

All MPI communication here uses MPI neighborhoods, meaning that it scales well with an increase of the number of used processes.

Thanks for the answer, I think that I got the info that I needed :smiling_face:. Since the coarse function is evaluated at the point then we are doing linear interpolation from the nearest nodes. My question is if it’s possible to consider more nodes in the coarse mesh and do a higher order interpolation? I hope I was clear now :sweat_smile: sorry for the confusion

Why would you need more nodes in the coarse grid?

The interpolation operator evaluates at the interpolation points of the finer grid, getting as good accuracy as it can get.

It won’t improve the accuracy of the interpolation by adding more nodes on the coarse grid, as it is still a linear function.

yes I understand what you mean, indeed what I would like to do is not usually done in FE, but we sometimes see that in FD. For instance, in 1D if we have three points x1<x2<x3 on the coarse grid we can do linear interpolation and estimate the function value at any point x in [x1,x2] by using only the data at x1 and x2. But if we use also the data at x3 we can estimate the function’s second derivative and “improve” the solution. For instance, 2nd order polynomials can be evaluated exactly :slight_smile:
But, as I said, I don’t know how we could do that easily with FE, since the grids are usually unstructured.

There is a fundamental difference between FD and FEM. FEM functions are defined at every point. If you use a linear basis, you can evaluate it at any point. The derivative is described by the derivatives of the basis functions multiplied by its coefficients. This means that you can have more accurate solutions by using higher order basis functions.

Hi all,
sorry for the late reply, I had to investigate a bit more… Indeed, what I would need in the FEM case is the restriction operator from the geometric multigrid. Is there any example where it is implemented? I was looking in the forum but didnt find anything
Thanks :slight_smile:

See for instance: Interpolation matrix with non matching meshes - #15 by bay_swiss

In general, the prolongation/restriction operator can be implicitly enforced using the non-matching interpolation.

Hi dokken,
thanks for the answer! I already tried with non-matching interpolation, but I wanted two try as well the restriction operator defined as the adjoint of the interpolation from coarse to fine. Is there a way for doing so?

Or, alternatively, a way of doing an orthogonal projection from fine to coarse? The result depends on the quadrature rule employed and should also be different than simple interpolation…

I tried to do so but I would need to use functions defined on different meshes in the same ufl form and I get the error TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'

A MWE is:

from mpi4py import MPI
from dolfinx import mesh, fem, io
import ufl
import numpy as np


def g(x):
    return 1.0 + np.sin(x[0]) ** 2 + 2 * x[1] ** 2


domain_fine = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.triangle)
domain_coarse = mesh.create_unit_square(MPI.COMM_WORLD, 16, 16, mesh.CellType.triangle)
V_fine = fem.FunctionSpace(domain_fine, ("CG", 1))
V_coarse = fem.FunctionSpace(domain_coarse, ("CG", 1))

# define some function g_f on the fine mesh
g_f = fem.Function(V_fine)
g_f.interpolate(g)

# try to do an ortogonal projection of g_f onto the coarse space
u = ufl.TrialFunction(V_coarse)
v = ufl.TestFunction(V_coarse)
m = u * v * ufl.dx
rhs = g_f * v * ufl.dx  # Error here: TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'
problem = fem.petsc.LinearProblem(m, rhs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
g_c_proj = problem.solve()

I there a way to overcome this?
Thanks :pray:

g_f is defined on the fine mesh while the problem itself is defined on the coarse mesh. You’d have to choose a scheme whereby g_f is interpolated or projected onto the coarse mesh accounting for the non matching meshes.

yes yes indeed what I would like to do is the L2 projection of g_f onto the coarse mesh

did you find how to do L2 projection of fine mesh solution to coarse mesh in fenicsx, I was wondering cause I am facing same trouble!!

Hey, unfortunately not. I’m considering using finite differences only. In any case, the focus of my research is time, not space.