Repository: fatiando/harmonica
Pull request #339 (2022-08-16 20:12:11)

Add a function to ignore the tesseroid with zero density or volume
When I am working on `tesseroid_layer` function, I had an error for tesseroid with zero volume. 
So I add a new checked in `tesseroid_gravity` that ignore the tesseroid with zero density or zero volume 


aguspesce: Hi @santisoler! 
I think that it is OK. Please, take a look at this! 

Patch:
diff --git a/harmonica/forward/_tesseroid_utils.py b/harmonica/forward/_tesseroid_utils.py
index f9e8914..eaa0b52 100644
--- a/harmonica/forward/_tesseroid_utils.py
+++ b/harmonica/forward/_tesseroid_utils.py
@@ -491,3 +491,37 @@ def _longitude_continuity(tesseroids):
     east[tess_to_be_changed] = ((east[tess_to_be_changed] + 180) % 360) - 180
     west[tess_to_be_changed] = ((west[tess_to_be_changed] + 180) % 360) - 180
     return tesseroids
+
+
+def _discard_null_tesseroids(tesseroids, density):
+    """
+    Discard tesseroid with zero volume or zero density
+
+    Parameters
+    ----------
+    tesseroids : 2d-array
+        Array containing the boundaries of the tesserois in the following
+        order: ``west``, ``east``, ``south``, ``north``, ``bottom``, ``top``
+        defined in a geocentric spherical coordinate system.
+    density : 1d-array
+        Array containing the density of each tesseroid in kg/m^3. Must have the
+        same size as the number of tesseroids.
+
+    Returns
+    -------
+    tesseroids : 2d-array
+        A copy of the ``tesseroids`` array that doesn't include the null
+        tesseroids (tesseroids with zero density or zero volume).
+    density : 1d-array
+        A copy of the ``density`` array that doesn't include the density values
+        for the null tesseroids (tesseroid with zero density or zero volume).
+    """
+    west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6))
+    # Mark prisms with zero volume as null prisms
+    null_tesseroids = (west == east) | (south == north) | (bottom == top)
+    # Mark prisms with zero density as null prisms
+    null_tesseroids[density == 0] = True
+    # Keep only non null prisms
+    tesseroids = tesseroids[np.logical_not(null_tesseroids), :]
+    density = density[np.logical_not(null_tesseroids)]
+    return tesseroids, density
diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py
index 3a3cbb4..02e8495 100644
--- a/harmonica/forward/tesseroid.py
+++ b/harmonica/forward/tesseroid.py
@@ -15,6 +15,7 @@ from ._tesseroid_utils import (
     _adaptive_discretization,
     _check_points_outside_tesseroids,
     _check_tesseroids,
+    _discard_null_tesseroids,
     gauss_legendre_quadrature,
     glq_nodes_weights,
 )
@@ -167,6 +168,8 @@ def tesseroid_gravity(
                 "Number of elements in density ({}) ".format(density.size)
                 + "mismatch the number of tesseroids ({})".format(tesseroids.shape[0])
             )
+        # Discard null tesseroids (zero density or zero volume)
+        tesseroids, density = _discard_null_tesseroids(tesseroids, density)
     # Get GLQ unscaled nodes, weights and number of nodes for each small
     # tesseroid
     glq_nodes, glq_weights = glq_nodes_weights(GLQ_DEGREES)


Test patch:
diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py
index 9db9533..feead05 100644
--- a/harmonica/tests/test_tesseroid.py
+++ b/harmonica/tests/test_tesseroid.py
@@ -15,6 +15,7 @@ from verde import grid_coordinates
 
 from ..constants import GRAVITATIONAL_CONST
 from ..forward._tesseroid_utils import (
+    _discard_null_tesseroids,
     _distance_tesseroid_point,
     _longitude_continuity,
     _split_tesseroid,
@@ -259,6 +260,41 @@ def test_stack_overflow_on_adaptive_discretization():
         )
 
 
+def test_discard_null_tesseroids():
+    """
+    Test if discarding invalid tesseroid works as expected
+    """
+    # Define a set of sample tesseroids including invalid ones
+    ellipsoid = boule.WGS84
+    top = ellipsoid.mean_radius
+    bottom = top - 1e3
+    tesseroids = np.array(
+        [
+            [-10, -5, -10, -5, bottom, top],  # ok tesseroid
+            [-10, -5, -5, 0, bottom, top],  # ok tesseroid (will set zero density)
+            [-10, -10, 0, 5, bottom, top],  # no volume due to easting boundaries
+            [-10, -5, 5, 5, bottom, top],  # no volume due to noting boundaries
+            [-5, 0, -10, -5, top, top],  # no volume due to radial boundaries
+            [-5, 0, -5, 0, bottom, top],  # ok tesseroid
+            [-5, -5, 5, 5, top, top],  # no volume due to multiple boundaries
+            [-5, 0, 5, 10, bottom, top],  # ok tesseroid
+        ]
+    )
+    densities = np.array([2400, 0, 2500, 2600, 2700, 2800, 2900, 3000])
+    tesseroids, densities = _discard_null_tesseroids(tesseroids, densities)
+    npt.assert_allclose(
+        tesseroids,
+        np.array(
+            [
+                [-10, -5, -10, -5, bottom, top],
+                [-5, 0, -5, 0, bottom, top],
+                [-5, 0, 5, 10, bottom, top],
+            ]
+        ),
+    )
+    npt.assert_allclose(densities, np.array([2400, 2800, 3000]))
+
+
 # --------------------------------------
 # Test tesseroid distance and dimensions
 # --------------------------------------