Cellwise interpolation of a discontuous function onto a (possibly nonmatching) mesh

dolfinx
1

Hey there,

I stumbled about the following problem last week and couldn’t resolve it on my own. I am working in a DG setting and want to interpolate a given function to a finer mesh. For that I am using a nonmatching interpolation, which works fine as long as the solution is continuous.

However, I found some behavior I could not explain, when I was trying to interpolate discontinuous functions. I then tried a cell-wise approach and the behavior looked even stranger. In the following MWE I tried to boil this down into a 1D case (later I’m more interested in a 2D version of that).

from mpi4py import MPI
import dolfinx as dfx
import numpy as np
import matplotlib.pyplot as plt

mesh_coarse = dfx.mesh.create_interval(MPI.COMM_WORLD, 4, (0, 1.0))
mesh_fine = dfx.mesh.create_interval(MPI.COMM_WORLD, 12, (0, 1.0))


fem_type = "DG"
V_h = dfx.fem.FunctionSpace(mesh_coarse, (fem_type, 1))
V_ref = dfx.fem.FunctionSpace(mesh_fine, (fem_type, 1))

u_h = dfx.fem.Function(V_h)
u_h.x.array[:] = np.arange(len(u_h.x.array))  # a coarse function with some jumps
V_h_x = np.array([coord[0] for coord in V_h.tabulate_dof_coordinates()])
V_ref_x = np.array([coord[0] for coord in V_ref.tabulate_dof_coordinates()])


u_h_int = dfx.fem.Function(V_ref)
u_h_glob_int = dfx.fem.Function(V_ref)
interpolation_data = dfx.fem.create_nonmatching_meshes_interpolation_data(
    V_ref.mesh._cpp_object, V_ref.element, V_h.mesh._cpp_object, padding=1e-9
)

u_h_glob_int.interpolate(u_h, nmm_interpolation_data=interpolation_data)

c_map = V_ref.mesh.topology.index_map(V_ref.mesh.topology.dim)
num_cells = c_map.size_local + c_map.num_ghosts
for cell in range(num_cells):
    u_h_int.interpolate(
        u_h,
        nmm_interpolation_data=interpolation_data,
        cells=np.array([cell]).astype(np.int32),
    )


plt.plot(V_h_x, u_h.x.array, marker="x", label="uh")
plt.plot(V_ref_x, u_h_int.x.array, marker="x", label="uh cellwise interpolated")
plt.plot(V_ref_x, u_h_glob_int.x.array, marker="x", label="uh interpolated")
plt.legend()
plt.show()

The code generates the following plot:
interpolation

I somehow want the green line to capture the jumps of the blue line, which is indeed possible in this setting (the fine mesh is a refinement of the coarse one). Note, that I plotted over coords[0], so the green line has the same amount of points as the orange one (my try on a cell-wise interpolation). As you see, the cell-wise approach didn’t work at all, and I’m not completely sure why. I also tried doing the interpolation a seperate Function instance and adding the relevant data to u_h_int, but the result was basically the same.

My two main questions are:

  • Is it possible to calculate a cell-wise interpolation in a DG context?
  • How can I capture jumps, when interpolating bewteen meshes?

Thanks a lot in advance :slight_smile:

2

DG-1 will be very tricky to capture accurately, as you will have multiple dofs at the same physical position, which the algorithm cannot logical distinguish.

If you want something discontinuous transfered, i would either:

  1. use a projection, as described here:
    L2 projection of fine mesh solution to a coarse mesh - #6 by dokken
  2. use DG-0, whose coordinate is at the cell center, and does not have duplicates.
3

Thanks a lot for the quick response. Would such a projection approach also work the other way around, i.e. projecting from a coarse to a finer mesh?

When calculating an error of some approximation against a reference solution this would feel more natural than the other way around.

4

Yes, that should work. Note that alignment between cells of the coarser and finer grid is crucial to avoid Gibbs as shown in Projection and interpolation — FEniCS Tutorial @ Sorbonne

1 Like