[jcapriot]

Summary

Fixes shift_to_discrete_topography putting locations off of the discrete surface on tree meshes.

PR Checklist

• If this is a work in progress PR, set as a Draft PR
• Linted my code according to the style guides.
• Added tests to verify changes to the code.
• Added necessary documentation to any new functions/classes following the
  expect style.
• Marked as ready for review (if this is was a draft PR), and converted
  to a Pull Request
• Tagged @simpeg/simpeg-developers when ready for review.

Reference issue

Came up while working on receiver locations for #1839, and follows up on the closest center vs containing cell comment in #1683 (#1683 (comment)).

What does this implement/fix?

There were a few ways a shifted location could end up off of the discrete surface on tree meshes, and on tensor meshes with uneven cell widths (like padding):

• With shift_horizontal=False, the elevation came from the surface cell with the closest center, which can be in the neighboring column (near a change in cell size where the surface steps, including the padding of a tensor mesh). It now walks down each location's own column to its highest active cell.
• With shift_horizontal=True and locations that include an elevation (which is what drape_electrodes_on_topography passes for 2D surveys), the horizontal position came from whatever cell held the original point, while the elevation came from the closest surface cell. Both now come from the closest surface cell, which is what horizontal only locations (and the old drapeTopotoLoc) already did.
• get_discrete_topography returned every active cell with an inactive cell directly above it, which included the floors of inactive pockets (and the ground under an overhang), so the horizontal shift could put a location underground. It now walks down through the inactive cells from the top of the mesh using each cell's neighbors, so it only finds the upper surface. Where finer inactive cells sit on a coarser surface cell, that cell's top is also split up into the horizontal positions of the finer cells above it (what the surface looks like from above).

For example, with the surface stepping down from z = 0 to z = -2 at x = 0:

import numpy as np
import discretize
from simpeg.utils import shift_to_discrete_topography

# surface at z = 0 (1 m cells) for x < 0, and z = -2 (2 m cells) for x > 0
mesh = discretize.TreeMesh([np.ones(32), np.ones(32)], "CC", diagonal_balance=True)
mesh.refine(
    lambda cell: mesh.max_level - (cell.center[0] > 0 or abs(cell.center[1]) > 2),
    diagonal_balance=True,
)
x, z = mesh.cell_centers.T
active_cells = z < np.where(x < 0, 0.0, -2.0)
electrode = np.array([[0.2, 0.0]])

print(shift_to_discrete_topography(mesh, electrode, active_cells, shift_horizontal=False))
print(shift_to_discrete_topography(mesh, electrode, active_cells, shift_horizontal=True))

On main this gives [0.2, 0.0] and [1.0, 0.0], which are both in the air (the surface is at z = -2 there). With this PR it gives [0.2, -2.0] and [-0.5, 0.0].

Shifted electrode locations on the same mesh, for a few more electrodes. Open circles are the original locations and filled circles the shifted ones, and the ringed ones end up in the air on main. Everything else lands in the same place.

The same thing on a tensor mesh, with a 10 m column (0 < x < 10) next to 1 m columns. The ringed locations end up in the ground on main, because the closest column center isn't the one they're in.

Both functions also accept integer active_cells now, and raise a clear error for an invalid topo_cell_cutoff on tree meshes (and, without the horizontal shift, for a column with no active cells).

Additional information

I'm not really sold on shift_horizontal=True as the default (especially for the nodal DC simulations), but I'll open a separate issue for that.

[codecov]

Codecov Report

❌ Patch coverage is 97.54601% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 86.85%. Comparing base (c593853) to head (1f7ef02).

Files with missing lines	Patch %	Lines
simpeg/utils/mesh_utils.py	94.20%	2 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1840      +/-   ##
==========================================
+ Coverage   86.81%   86.85%   +0.03%     
==========================================
  Files         436      436              
  Lines       57838    57952     +114     
  Branches     5457     5460       +3     
==========================================
+ Hits        50214    50332     +118     
+ Misses       6187     6185       -2     
+ Partials     1437     1435       -2     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.