diff --git a/simpeg_drivers/components/topography.py b/simpeg_drivers/components/topography.py index f94d351c..5030b7e1 100644 --- a/simpeg_drivers/components/topography.py +++ b/simpeg_drivers/components/topography.py @@ -28,6 +28,7 @@ from geoh5py.data import NumericData from geoh5py.objects.surveys.electromagnetics.base import LargeLoopGroundEMSurvey from geoh5py.shared import Entity +from scipy.spatial import cKDTree from simpeg_drivers.components.data import InversionData from simpeg_drivers.components.locations import InversionLocations @@ -37,7 +38,6 @@ active_from_xyz, floating_active, get_containing_cells, - get_neighbouring_cells, mask_vertices_and_cells, ) @@ -121,7 +121,7 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray: active_cells = active_from_xyz( mesh.entity, vertices, - grid_reference="bottom" if forced_to_surface else "center", + grid_reference="center", triangulation=cells, ) @@ -181,22 +181,14 @@ def expand_actives( containing_cells = get_containing_cells(mesh.mesh, data) active_cells[containing_cells] = True - # Apply extra active cells to ensure connectivity for tree meshes - if isinstance(mesh.mesh, TreeMesh): - neighbours = get_neighbouring_cells(mesh.mesh, containing_cells) - neighbours_xy = np.r_[neighbours[0] + neighbours[1]] - - neighbours_xy = neighbours_xy[neighbours_xy != -1] - # Make sure the new actives are connected to the old actives - new_actives = ~active_cells[neighbours_xy] - if np.any(new_actives): - neighbours = get_neighbouring_cells( - mesh.mesh, neighbours_xy[new_actives] - ) - neighbours_z = np.r_[neighbours[2][0]] - neighbours_z = neighbours_z[neighbours_z != -1] - active_cells[neighbours_z] = True # z-axis neighbours - - active_cells[neighbours_xy] = True # xy-axis neighbours + # Apply extra active cells to ensure connectivity for neighbours + tree = cKDTree(mesh.mesh.cell_centers[containing_cells]) + rad, ind = tree.query(mesh.mesh.cell_centers) + neighbours_xy = rad < (3 * mesh.mesh.h[0].min()) + neighbours_xy &= ( + mesh.mesh.cell_centers[containing_cells, :][ind, -1] + >= mesh.mesh.cell_centers[:, -1] + ) + active_cells[neighbours_xy] = True # xy-axis neighbours return active_cells diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index d26e324c..7100f55e 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -505,7 +505,11 @@ def get_containing_cells( inds = np.r_[inds, np.hstack(line_ind)] elif isinstance(mesh, TensorMesh): - locations = data.drape_locations(np.unique(data.locations, axis=0)) + potentials = data.entity.vertices + currents = data.entity.current_electrodes.vertices + locations = np.unique(np.r_[potentials, currents], axis=0) + + locations = data.drape_locations(np.unique(locations, axis=0)) xi = np.searchsorted(mesh.nodes_x, locations[:, 0]) - 1 yi = np.searchsorted(mesh.nodes_y, locations[:, -1]) - 1 inds = xi + yi * mesh.shape_cells[0] diff --git a/tests/run_tests/driver_dc_2d_test.py b/tests/run_tests/driver_dc_2d_test.py index bc2d7d78..5af1d0fc 100644 --- a/tests/run_tests/driver_dc_2d_test.py +++ b/tests/run_tests/driver_dc_2d_test.py @@ -43,7 +43,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 11.136935742296085, "phi_d": 5570, "phi_m": 314} +target_run = {"data_norm": 10.945383968572745, "phi_d": 5210, "phi_m": 340} def test_dc_2d_fwr_run( diff --git a/tests/run_tests/driver_dc_test.py b/tests/run_tests/driver_dc_test.py index 9e40ca8d..0f42327d 100644 --- a/tests/run_tests/driver_dc_test.py +++ b/tests/run_tests/driver_dc_test.py @@ -38,7 +38,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 0.15320935486917722, "phi_d": 25.7, "phi_m": 3580} +target_run = {"data_norm": 0.14272756694409652, "phi_d": 12.6, "phi_m": 3580} def test_dc_3d_fwr_run( diff --git a/tests/run_tests/driver_ip_2d_test.py b/tests/run_tests/driver_ip_2d_test.py index 02394521..2ea4b34b 100644 --- a/tests/run_tests/driver_ip_2d_test.py +++ b/tests/run_tests/driver_ip_2d_test.py @@ -39,7 +39,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 0.1244717397585979, "phi_d": 15500, "phi_m": 0.0002845} +target_run = {"data_norm": 0.12345083692791135, "phi_d": 15700, "phi_m": 0.000206} def test_ip_2d_fwr_run( diff --git a/tests/run_tests/driver_ip_test.py b/tests/run_tests/driver_ip_test.py index 5a292507..0d019035 100644 --- a/tests/run_tests/driver_ip_test.py +++ b/tests/run_tests/driver_ip_test.py @@ -37,7 +37,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 0.005219568872697125, "phi_d": 680, "phi_m": 1.01e-6} +target_run = {"data_norm": 0.004747013060439004, "phi_d": 562, "phi_m": 1.43e-06} def test_ip_3d_fwr_run( diff --git a/tests/run_tests/driver_joint_cross_gradient_test.py b/tests/run_tests/driver_joint_cross_gradient_test.py index ae3e815f..d014cd23 100644 --- a/tests/run_tests/driver_joint_cross_gradient_test.py +++ b/tests/run_tests/driver_joint_cross_gradient_test.py @@ -325,13 +325,13 @@ def test_joint_cross_gradient_inv_run( # the scaling from its total misfit. np.testing.assert_allclose( driver.directives.scale_misfits.scalings, - [1, 0.7558, 0.5, 0.5, 0.6710], + [1.0, 0.755269, 0.5, 0.5, 0.675174], atol=1e-3, ) # Check that scaling * chi factor is reflected in data misfit multipliers np.testing.assert_allclose( driver.data_misfit.multipliers, - [0.8, 0.6046, 0.5, 0.5, 0.6710], + [0.8, 0.604215, 0.5, 0.5, 0.675174], atol=1e-3, ) diff --git a/tests/run_tests/driver_mt_test.py b/tests/run_tests/driver_mt_test.py index b0a68eac..c49afc7d 100644 --- a/tests/run_tests/driver_mt_test.py +++ b/tests/run_tests/driver_mt_test.py @@ -42,7 +42,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 0.025087759238073448, "phi_d": 0.385, "phi_m": 3} +target_run = {"data_norm": 0.003584600661140228, "phi_d": 4.45, "phi_m": 5.56} def setup_data(workspace, survey):