Interpolating scalar/vector fields between fine and coarse meshes in GetFEM++

[Djoumessi]

Hello everyone,

I’m working on a multiphysics problem where I’m coupling a reaction-diffusion system (which requires a fine mesh) with a mechanical model (which typically converges well on a coarser mesh). To reduce computational cost, I want to solve each system on a separate mesh and then transfer fields between the meshes, especially from the fine one to the coarse one and vice versa.

I’ve been using gf.compute_interpolate_on and gf_compute_extrapolate_on to transfer scalar and vector fields (e.g., concentration, internal strain, etc.) between the meshes. However, I’m facing difficulties in getting this to work reliably when the source and target meshes are different (though they cover the same domain), especially:

• When interpolating from a fine mesh (high resolution) to a coarser one for the mechanical model;
• When extrapolating from a coarse mechanical mesh to a fine reaction-diffusion mesh to feed back deformation or strain information.

In some tests, the interpolation seemed to return NaNs or values at zero, likely due to the lack of correspondence in mesh nodes or insufficient coverage. I suspect this is because the meshes are not nested and the interpolation might miss points outside the convex hull.

Here’s a simplified version of what I’m trying to do:

phi_fine = md_diffusion.variable("phi")  # solved on fine mesh
phi_coarse = gf.compute_interpolate_on(mf_fine, phi_fine, mf_coarse)  # interpolation
md_mechanics.add_fem_data("phi", mf_coarse)
md_mechanics.set_variable("phi", phi_coarse)

Similarly, I tried using extrapolation from coarse to fine, but encountered instability or inaccuracies.

I’m aware that interpolating between different FEM spaces (e.g., P2 to P1 or P1 to P2) on the same mesh works well, but I need to extend this to different meshes for efficiency.

My main goal is to make this kind of coupling efficient:

• fine mesh for the physics that needs it (e.g., diffusion),
• coarse mesh for what can afford it (e.g., mechanics),
• and reliable, accurate transfer between them.

Has anyone encountered a similar issue or implemented a robust solution for interpolation/extrapolation across non-nested meshes in GetFEM++?

Thanks in advance for your help and suggestions!

Best regards,

[Konstantinos.Poulios]

Hi,
compute_interpolate_on is one of the legacy functions. These days we use model.interpolation() instead.

What might be more interesting for you is to couple the equations between the two meshes instead of interpolating on the fly. For example, I use the following syntax to couple two different physics defined on a 1d mesh m1 and a 2d mesh m2:

md.add_interpolate_transformation_from_expression("surface", m2, m1, "X")
md.add_interpolate_transformation_from_expression("substrate", m1, m2, "X")
md.add_macro("H", "H0+sqr(X(1))/2-Interpolate(U,substrate)(2)")

Another useful operation is projecting between 2 FE spaces. You can study and learn something useful from the following snippet:

mim0 = gf.MeshIm(m0, 3)
X0 = np.load("x_1921.npy")
if extra_refine:
  mfx0c = gf.MeshFem(m0c, 1)
  mfx0c.set_classical_fem(1)
  mdaux0 = gf.Model("real")
  mdaux0.add_fem_variable("x", mfx0)
  mdaux0.add_initialized_fem_data("xc", mfx0c, X0)
  mdaux0.add_interpolate_transformation_from_expression("coarse2fine", m0, m0c, "X")
  X = gf.linsolve_mumps(gf.asm_mass_matrix(mim0, mfx0),
                        gf.asm_generic(mim0, 1, "Interpolate(xc,coarse2fine)*Test_x", -1, mdaux0)).flatten()
else:
  X = X0
mfx0.export_to_vtu("levelset.vtu", mfx0, X, "x")

[Djoumessi]

Dear Konstantinos,

Thank you very much for your kind help.

I ended up implementing the second approach, which is an 𝐿2 projection, consists in solving the following equation or something similar:
\int_{\Omega} \phi_{\text{fine}} \cdot v_{\text{coarse}} \, dx = \int_{\Omega} \tilde{\phi}_{\text{coarse}} \cdot v_{\text{coarse}} \, dx
The method works very well overall. However, I noticed that the following line:

md.add_interpolate_transformation_from_expression("source2target", m_target, m_source, "X")

only works reliably when the coarse mesh has a sufficiently dense distribution of nodes relative to the fine mesh. Otherwise, it seems that the transformation fails, possibly because some nodes on the target mesh cannot be mapped to any element of the source mesh.

Please correct me if I’m wrong.

Thanks again for your support.

Best regards,

[Konstantinos.Poulios]

Can you show a region of the two overlapping meshes to get an idea?

If you use the interpolation in an L2 projection, then interpolation will never occur on nodes, only on integration points.

GWFL has a syntax for dealing with failed interpolation, read about interpolation_filter in:

Compute arbitrary terms - high-level generic assembly procedures - Generic...

Not 100% sure but I guess the interpolate_transformation_from_expression must be compatible with filters. If you use the filter you need to decide yourself what to do with failed points. One option is to just ignore them if they do not matter.

[Djoumessi]

Dear Konstantinos,

Thank you very much for your helpful and detailed response.

Indeed, I chose to use an L² projection because I expected it to be a more stable and general approach, even though the interpolation is performed at the quadrature points rather than directly at the mesh nodes — as you correctly pointed out.

Regarding the two meshes (shown in the attached image), they actually represent the same computational domain, but with different mesh sizes. The fine mesh captures more geometric detail, while the coarse one is used for efficiency or multiscale coupling.

I haven’t tested the use of Interpolate_filter yet, but I will explore it and let you know how it performs in practice. It seems like a promising solution to handle failed interpolation cases more gracefully.

Thanks again for your support and guidance!

Best regards,