BUG: Incorrect logic for zero_outside=True interpolation on TreeMesh for cell_centered, edge, and face interpolation matrices.

[dthillaithevan]

Describe the issue:

maps.Mesh2Mesh appears to be giving incorrect results when interpolating from a TreeMesh onto another TreeMesh.

In the example below I am attempting to interpolate the cell volumes and a homogeneous model of 1's onto a fine mesh. In both cases the interpolated models on the fine mesh is giving values that appear to be incorrect. E.g. In the homogeneous case I would expect a resulting interpolated model of ones but I am getting a constant model filled with 8.

Reproducible code example:

import discretize
from simpeg import maps, utils
import matplotlib.pyplot as plt
from simpeg.utils import mkvc
from discretize import TreeMesh

import numpy as np 

def create_tree_mesh(dh: float):
    # Modified from https://discretize.simpeg.xyz/en/main/tutorials/mesh_generation/4_tree_mesh.html
    dx = dy = dz = dh

    x_length = 300.0  # domain width in x
    y_length = 300.0  # domain width in y
    z_length = 300.0  # domain width in y
    
    # Compute number of base mesh cells required in x and y
    nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0)))
    nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0)))
    nbcz = 2 ** int(np.round(np.log(z_length / dz) / np.log(2.0)))
    
    # Define the base mesh
    hx = [(dx, nbcx)]
    hy = [(dy, nbcy)]
    hz = [(dz, nbcz)]
    mesh = TreeMesh([hx, hy, hz], x0="CCC", diagonal_balance=True)
    
    # Refine surface topography
    [xx, yy] = np.meshgrid(mesh.nodes_x, mesh.nodes_y)
    zz = -3 * np.exp((xx**2 + yy**2) / 100**2) + 50.0
    pts = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)]
    padding = [[0, 0, 2], [0, 0, 2]]
    mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False)
    
    # Refine box
    xp, yp, zp = np.meshgrid([-40.0, 40.0], [-40.0, 40.0], [-60.0, 0.0])
    xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
    mesh.refine_bounding_box(xyz, padding_cells_by_level=padding, finalize=False)
    mesh.finalize()

    return mesh

def plot_slice(v: np.ndarray, mesh: TreeMesh):
    fig = plt.figure(figsize=(6, 6))
    ax = fig.add_subplot(111)
    
    if mesh.dim == 3:
        mesh.plot_slice(np.log10(v), normal="Y", ax=ax, ind=int(mesh.h[1].size / 2), grid=True)
    elif mesh.dim == 2:
        mesh.plot_image(np.log10(v), grid=True, ax = ax)
    ax.set_title("Cell Log-Volumes at Y = 0 m")
    
    plt.show()

if __name__ == "__main__":
    
    # Coarse mesh
    mesh_c = create_tree_mesh(10)
    # Finer mesh
    mesh_f = create_tree_mesh(5)
    
    # Coarse mesh cell volumnes
    coarse_vols = mesh_c.cell_volumes
    coarse_ones = np.ones_like(coarse_vols)
    
    # Mapping from coarse to fine mesh
    mod = maps.Mesh2Mesh([mesh_f, mesh_c])
    # Interpolate coarse volumes onto fine mesh
    mapped_coarse_vols = mod * coarse_vols
    
    # Interpolate homogeneous field (=1) onto fine mesh
    mapped_coarse_ones = mod * coarse_ones
    
    plot_slice(coarse_vols, mesh_c)    
    plot_slice(mapped_coarse_vols, mesh_f)
    print (coarse_vols.max(), coarse_vols.min()) # 64000.0, 1000.0
    print (mapped_coarse_vols.max(), mapped_coarse_vols.min()) # 8000.0 8000.0
    
    
    plot_slice(coarse_ones, mesh_c)    
    plot_slice(mapped_coarse_ones, mesh_f)
    print (coarse_ones.max(), coarse_ones.min()) # 1.0, 1.0
    print (mapped_coarse_ones.max(), mapped_coarse_ones.min()) # 8.0, 8.0

Error message:

Runtime information:

--------------------------------------------------------------------------------
  Date: Tue Nov 18 09:30:51 2025 UTC

                OS : Darwin (macOS 15.5)
            CPU(s) : 12
           Machine : arm64
      Architecture : 64bit
               RAM : 48.0 GiB
       Environment : IPython
       File system : apfs

  Python 3.12.12 | packaged by conda-forge | (main, Oct 13 2025, 14:47:45)
  [Clang 19.1.7 ]

            simpeg : 0.24.0
        discretize : 0.12.0
       pymatsolver : 0.3.1
             numpy : 2.3.3
             scipy : 1.16.2
        matplotlib : 3.10.7
            geoana : 0.8.0
            libdlf : 0.3.0
             numba : 0.62.1
           sklearn : 1.7.2
            pandas : 2.3.3
           IPython : 9.6.0
            plotly : 6.3.1
               vtk : 9.5.2
            choclo : 0.3.2
--------------------------------------------------------------------------------

Context for the issue:

No response

[santisoler]

Hi @dthillaithevan. Thanks for opening this issue. I think there's a bug indeed. I don't think it's coming directly form the maps.Mesh2Mesh directive, but from the TreeMesh.get_interpolation_matrix method.

I took your example and instead of using the mapping, I'm manually creating the interpolation matrix and multiply it by the cell volumes of the coarser mesh:

from simpeg.utils import mkvc
from discretize import TreeMesh

import numpy as np


def create_tree_mesh(dh: float):
    # Modified from https://discretize.simpeg.xyz/en/main/tutorials/mesh_generation/4_tree_mesh.html
    dx = dy = dz = dh

    x_length = 300.0  # domain width in x
    y_length = 300.0  # domain width in y
    z_length = 300.0  # domain width in y

    # Compute number of base mesh cells required in x and y
    nbcx = 2 ** int(np.round(np.log(x_length / dx) / np.log(2.0)))
    nbcy = 2 ** int(np.round(np.log(y_length / dy) / np.log(2.0)))
    nbcz = 2 ** int(np.round(np.log(z_length / dz) / np.log(2.0)))

    # Define the base mesh
    hx = [(dx, nbcx)]
    hy = [(dy, nbcy)]
    hz = [(dz, nbcz)]
    mesh = TreeMesh([hx, hy, hz], x0="CCC", diagonal_balance=True)

    # Refine surface topography
    [xx, yy] = np.meshgrid(mesh.nodes_x, mesh.nodes_y)
    zz = -3 * np.exp((xx**2 + yy**2) / 100**2) + 50.0
    pts = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)]
    padding = [[0, 0, 2], [0, 0, 2]]
    mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False)

    # Refine box
    xp, yp, zp = np.meshgrid([-40.0, 40.0], [-40.0, 40.0], [-60.0, 0.0])
    xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)]
    mesh.refine_bounding_box(xyz, padding_cells_by_level=padding, finalize=False)
    mesh.finalize()

    return mesh


if __name__ == "__main__":

    mesh_coarse = create_tree_mesh(10)
    mesh_fine = create_tree_mesh(5)

    # Build interpolation matrix
    interpolation_matrix = mesh_coarse.get_interpolation_matrix(
        mesh_fine.cell_centers,
        location_type="cell_centers",
        zeros_outside=True,
    )

    result = interpolation_matrix @ mesh_coarse.cell_volumes

    print("Coarse mesh cell volumes")
    print(mesh_coarse.cell_volumes.min(), mesh_coarse.cell_volumes.max())

    print("Interpolated cell volumes")
    print(result.min(), result.max())

Coarse mesh cell volumes
1000.0 64000.0
Interpolated cell volumes
8000.0 8000.0

So I think the issue is actually in discretize, and not in SimPEG. I suspect the interpolation matrix might be including cell volumes, but I could be wrong...

@jcapriot, do you have a better idea of what's going on here?

[santisoler]

As a side note, just noted that the TensorMesh.get_interpolation_matrix takes loc as first argument, while TreeMesh.get_interpolation_matrix takes locs instead. Minor difference, but can bite us when passing arguments with keys.

https://discretize.simpeg.xyz/en/main/api/generated/discretize.TensorMesh.get_interpolation_matrix.html#discretize.TensorMesh.get_interpolation_matrix

https://discretize.simpeg.xyz/en/main/api/generated/discretize.TreeMesh.get_interpolation_matrix.html#discretize.TreeMesh.get_interpolation_matrix

[jcapriot]

Interesting bug! We don't use the zero_outside=True option very often. Looks like this bug got introduced in #277, 3 years ago!

[dthillaithevan]

Thanks @santisoler and @jcapriot for looking into this!

I had suspected it was something to do with cell volumes being propagated incorrectly or not being normalised, due to the factor of 8 showing up across different interpolations, but I am afraid I am not smart enough to debug this on my own, so I appreciate your help on this.

For the time being, I have written a function which manually approximates the interpolation matrix which is doing the trick for now :)