Merge pull request #838 from bp/none_faces_issues

further work on find faces for curtains

resqpy/grid_surface/_find_faces.py

@@ -1145,8 +1145,12 @@ def find_faces_to_represent_surface_regular_optimised(grid,
     if progress_fn is not None:
         progress_fn(0.9)
 
-    assert k_faces_kji0 is not None or j_faces_kji0 is not None or i_faces_kji0 is not None,  \
-        f'did not find any faces to represent {name}: surface does not intersect grid?'
+    if k_faces_kji0 is None and j_faces_kji0 is None and i_faces_kji0 is None:
+        log.error(f'did not find any faces to represent {name}: surface does not intersect grid?')
+        if return_properties:
+            return (None, {})
+        else:
+            return None
 
     log.debug("converting face sets into grid connection set")
     # NB: kji0 arrays in internal face protocol: used as cell_kji0 with polarity of 1
@@ -1210,14 +1214,26 @@ def find_faces_to_represent_surface_regular_optimised(grid,
     if return_bisector:
         if is_curtain:
             log.debug("preparing columns bisector")
-            bisector = column_bisector_from_face_indices((grid.nj, grid.ni), j_faces_kji0[:, 1:], i_faces_kji0[:, 1:])
+            if j_faces_kji0 is None:
+                j_faces_ji0 = np.empty((0, 2), dtype = np.int32)
+            else:
+                j_faces_ji0 = j_faces_kji0[:, 1:]
+            if i_faces_kji0 is None:
+                i_faces_ji0 = np.empty((0, 2), dtype = np.int32)
+            else:
+                i_faces_ji0 = i_faces_kji0[:, 1:]
+            bisector = column_bisector_from_face_indices((grid.nj, grid.ni), j_faces_ji0, i_faces_ji0)
             # log.debug('finished preparing columns bisector')
         else:
             log.debug("preparing cells bisector")
-            bisector, is_curtain = bisector_from_face_indices(tuple(grid.extent_kji), k_faces_kji0, j_faces_kji0,
-                                                              i_faces_kji0, raw_bisector)
-            if is_curtain:
-                bisector = bisector[0]  # reduce to a columns property
+            if k_faces_kji0 is None and j_faces_kji0 is None and i_faces_kji0 is None:
+                bisector = np.ones((grid.nj, grid.ni), dtype = bool)
+                is_curtain = True
+            else:
+                bisector, is_curtain = bisector_from_face_indices(tuple(grid.extent_kji), k_faces_kji0, j_faces_kji0,
+                                                                  i_faces_kji0, raw_bisector)
+                if is_curtain:
+                    bisector = bisector[0]  # reduce to a columns property
 
     # note: following is a grid cells property, not a gcs property
     if return_shadow:

tests/unit_tests/conftest.py

@@ -533,8 +533,63 @@ def small_grid_and_surface_no_k(tmp_model: Model) -> Tuple[grr.RegularGrid, rqs.
     grid = grr.RegularGrid(tmp_model, extent_kji = extent_kji, dxyz = dxyz, crs_uuid = crs_uuid, title = title)
     grid.create_xml()
 
-    n_points = 100
-    points = np.array(([0.37, 0.043, 0.017], [0.37, 0.043, 0.921], [0.61, 0.903, 0.027], [0.61, 0.903, 0.981]),
+    points = np.array(([0.37, -0.043, -0.017], [0.37, -0.043, 1.021], [0.61, 1.003, -0.027], [0.61, 1.003, 1.081]),
+                      dtype = float) * extent
+    triangles = np.array(([0, 1, 2], [1, 2, 3]), dtype = np.int32)
+    surface = rqs.Surface(tmp_model, crs_uuid = crs_uuid, title = "small_curtain")
+    surface.set_from_triangles_and_points(triangles, points)
+    surface.triangles_and_points()
+    surface.write_hdf5()
+    surface.create_xml()
+
+    tmp_model.store_epc()
+
+    return grid, surface
+
+
+@pytest.fixture
+def small_grid_and_i_curtain_surface(tmp_model: Model) -> Tuple[grr.RegularGrid, rqs.Surface]:
+    """Creates a small RegularGrid and a curtain triangular surface."""
+    crs = Crs(tmp_model)
+    crs.create_xml()
+
+    extent = 10
+    extent_kji = (extent, extent, extent)
+    dxyz = (1.0, 1.0, 1.0)
+    crs_uuid = crs.uuid
+    title = "small_grid"
+    grid = grr.RegularGrid(tmp_model, extent_kji = extent_kji, dxyz = dxyz, crs_uuid = crs_uuid, title = title)
+    grid.create_xml()
+
+    points = np.array(([0.47, -0.043, -0.017], [0.47, -0.043, 1.021], [0.47, 1.003, -0.027], [0.47, 1.003, 1.081]),
+                      dtype = float) * extent
+    triangles = np.array(([0, 1, 2], [1, 2, 3]), dtype = np.int32)
+    surface = rqs.Surface(tmp_model, crs_uuid = crs_uuid, title = "small_curtain")
+    surface.set_from_triangles_and_points(triangles, points)
+    surface.triangles_and_points()
+    surface.write_hdf5()
+    surface.create_xml()
+
+    tmp_model.store_epc()
+
+    return grid, surface
+
+
+@pytest.fixture
+def small_grid_and_j_curtain_surface(tmp_model: Model) -> Tuple[grr.RegularGrid, rqs.Surface]:
+    """Creates a small RegularGrid and a curtain triangular surface."""
+    crs = Crs(tmp_model)
+    crs.create_xml()
+
+    extent = 10
+    extent_kji = (extent, extent, extent)
+    dxyz = (1.0, 1.0, 1.0)
+    crs_uuid = crs.uuid
+    title = "small_grid"
+    grid = grr.RegularGrid(tmp_model, extent_kji = extent_kji, dxyz = dxyz, crs_uuid = crs_uuid, title = title)
+    grid.create_xml()
+
+    points = np.array(([-0.043, 0.59, -0.017], [-0.043, 0.59, 1.021], [1.003, 0.59, -0.027], [1.003, 0.59, 1.081]),
                       dtype = float) * extent
     triangles = np.array(([0, 1, 2], [1, 2, 3]), dtype = np.int32)
     surface = rqs.Surface(tmp_model, crs_uuid = crs_uuid, title = "small_curtain")
@@ -548,6 +603,33 @@ def small_grid_and_surface_no_k(tmp_model: Model) -> Tuple[grr.RegularGrid, rqs.
     return grid, surface
 
 
+@pytest.fixture
+def small_grid_and_missing_surface(tmp_model: Model) -> Tuple[grr.RegularGrid, rqs.Surface]:
+    """Creates a small RegularGrid and a triangular surface that does not intersect the grid."""
+    crs = Crs(tmp_model)
+    crs.create_xml()
+
+    extent = 10
+    extent_kji = (extent, extent, extent)
+    dxyz = (1.0, 1.0, 1.0)
+    crs_uuid = crs.uuid
+    title = "small_grid"
+    grid = grr.RegularGrid(tmp_model, extent_kji = extent_kji, dxyz = dxyz, crs_uuid = crs_uuid, title = title)
+    grid.create_xml()
+
+    points = np.array(([0.5, 0.0, 3.0], [0.0, 0.5, 3.0], [-3.0, 0.0, 0.0], [0.0, -3.0, 0.0]), dtype = float) * extent
+    triangles = np.array(([0, 1, 2], [1, 2, 3]), dtype = np.int32)
+    surface = rqs.Surface(tmp_model, crs_uuid = crs_uuid, title = "small_missing")
+    surface.set_from_triangles_and_points(triangles, points)
+    surface.triangles_and_points()
+    surface.write_hdf5()
+    surface.create_xml()
+
+    tmp_model.store_epc()
+
+    return grid, surface
+
+
 @pytest.fixture
 def small_grid_and_surface_no_j(tmp_model: Model) -> Tuple[grr.RegularGrid, rqs.Surface]:
     """Creates a small RegularGrid and tiny triangular surface which will map to no J faces."""

tests/unit_tests/test_grid_surface.py

@@ -159,6 +159,86 @@ def test_find_faces_to_represent_surface_no_i_faces(small_grid_and_surface_no_i)
     np.testing.assert_array_equal(fip_normal, fip_old_fuddy_duddy)
 
 
+def test_find_faces_to_represent_surface_no_i_or_k_faces(small_grid_and_j_curtain_surface):
+    # Arrange
+    grid = small_grid_and_j_curtain_surface[0]
+    surface = small_grid_and_j_curtain_surface[1]
+    old_fuddy_duddy_crs = rqc.Crs(surface.model, xy_units = 'ft', z_units = 'chain')
+    old_fuddy_duddy_crs.create_xml()
+    surface.model.store_epc()
+    s2 = rqs.Surface(surface.model, uuid = surface.uuid)
+    s2.change_crs(old_fuddy_duddy_crs)
+    name = "test"
+    assert grid.is_aligned
+
+    # Act
+    gcs_normal = rqgs.find_faces_to_represent_surface_regular(grid, surface, name)
+    cip_normal = gcs_normal.cell_index_pairs
+    fip_normal = gcs_normal.face_index_pairs
+
+    gcs_dense = rqgs.find_faces_to_represent_surface_regular_dense_optimised(grid, surface, name)
+    cip_dense = gcs_dense.cell_index_pairs
+    fip_dense = gcs_dense.face_index_pairs
+
+    gcs_optimised = rqgs.find_faces_to_represent_surface_regular_optimised(grid, surface, name)
+    cip_optimised = gcs_optimised.cell_index_pairs
+    fip_optimised = gcs_optimised.face_index_pairs
+
+    gcs_old_fuddy_duddy = rqgs.find_faces_to_represent_surface_regular_optimised(grid, s2, name)
+    cip_old_fuddy_duddy = gcs_old_fuddy_duddy.cell_index_pairs
+    fip_old_fuddy_duddy = gcs_old_fuddy_duddy.face_index_pairs
+
+    # Assert – quite harsh as gcs face ordering could legitimately vary
+    np.testing.assert_array_equal(cip_normal, cip_dense)
+    np.testing.assert_array_equal(fip_normal, fip_dense)
+    np.testing.assert_array_equal(cip_normal, cip_optimised)
+    np.testing.assert_array_equal(fip_normal, fip_optimised)
+    np.testing.assert_array_equal(cip_normal, cip_old_fuddy_duddy)
+    np.testing.assert_array_equal(fip_normal, fip_old_fuddy_duddy)
+    _, fip = gcs_optimised.list_of_cell_face_pairs_for_feature_index()
+    assert np.all(fip[:, :, 0] == 1)  # all J faces
+
+
+def test_find_faces_to_represent_surface_no_j_or_k_faces(small_grid_and_i_curtain_surface):
+    # Arrange
+    grid = small_grid_and_i_curtain_surface[0]
+    surface = small_grid_and_i_curtain_surface[1]
+    old_fuddy_duddy_crs = rqc.Crs(surface.model, xy_units = 'ft', z_units = 'chain')
+    old_fuddy_duddy_crs.create_xml()
+    surface.model.store_epc()
+    s2 = rqs.Surface(surface.model, uuid = surface.uuid)
+    s2.change_crs(old_fuddy_duddy_crs)
+    name = "test"
+    assert grid.is_aligned
+
+    # Act
+    gcs_normal = rqgs.find_faces_to_represent_surface_regular(grid, surface, name)
+    cip_normal = gcs_normal.cell_index_pairs
+    fip_normal = gcs_normal.face_index_pairs
+
+    gcs_dense = rqgs.find_faces_to_represent_surface_regular_dense_optimised(grid, surface, name)
+    cip_dense = gcs_dense.cell_index_pairs
+    fip_dense = gcs_dense.face_index_pairs
+
+    gcs_optimised = rqgs.find_faces_to_represent_surface_regular_optimised(grid, surface, name)
+    cip_optimised = gcs_optimised.cell_index_pairs
+    fip_optimised = gcs_optimised.face_index_pairs
+
+    gcs_old_fuddy_duddy = rqgs.find_faces_to_represent_surface_regular_optimised(grid, s2, name)
+    cip_old_fuddy_duddy = gcs_old_fuddy_duddy.cell_index_pairs
+    fip_old_fuddy_duddy = gcs_old_fuddy_duddy.face_index_pairs
+
+    # Assert – quite harsh as gcs face ordering could legitimately vary
+    np.testing.assert_array_equal(cip_normal, cip_dense)
+    np.testing.assert_array_equal(fip_normal, fip_dense)
+    np.testing.assert_array_equal(cip_normal, cip_optimised)
+    np.testing.assert_array_equal(fip_normal, fip_optimised)
+    np.testing.assert_array_equal(cip_normal, cip_old_fuddy_duddy)
+    np.testing.assert_array_equal(fip_normal, fip_old_fuddy_duddy)
+    _, fip = gcs_optimised.list_of_cell_face_pairs_for_feature_index()
+    assert np.all(fip[:, :, 0] == 2)  # all I faces
+
+
 def test_find_faces_to_represent_surface_regular_optimised_constant_agitation(small_grid_and_surface):
     # Arrange
     grid = small_grid_and_surface[0]
@@ -309,6 +389,7 @@ def test_find_faces_to_represent_curtain_regular_optimised_with_return_propertie
     return_properties.append("depth")
     return_properties.append("triangle")
     return_properties.append("flange bool")
+    return_properties.append("grid bisector")
     (
         gcs_optimised,
         properties_optimised,
@@ -322,6 +403,7 @@ def test_find_faces_to_represent_curtain_regular_optimised_with_return_propertie
     depths_optimised = properties_optimised["depth"]
     offsets_optimised = properties_optimised["offset"]
     flange_optimised = properties_optimised["flange bool"]
+    bisector_optimised, is_curtain_optimised = properties_optimised["grid bisector"]
 
     # Assert – quite harsh as faces could legitimately be in different order
     np.testing.assert_array_equal(cip_normal, cip_optimised)
@@ -335,6 +417,8 @@ def test_find_faces_to_represent_curtain_regular_optimised_with_return_propertie
     assert np.all(triangles_optimised >= 0)
     assert flange_optimised.shape == offsets_optimised.shape
     assert not np.any(flange_optimised)
+    assert bisector_optimised.shape == (grid.nj, grid.ni)
+    assert is_curtain_optimised
 
 
 def test_find_faces_to_represent_surface_regular_dense_optimised_with_return_properties(small_grid_and_surface,):
@@ -640,3 +724,22 @@ def test_shadow_from_faces_flat_surface_k_hole():
         [[2, 2, 2], [2, 0, 2], [3, 2, 2]],
         [[2, 2, 2], [2, 0, 2], [2, 2, 2]],
     ]))
+
+
+def test_find_faces_to_represent_surface_missing_grid(small_grid_and_missing_surface):
+    # Arrange
+    grid = small_grid_and_missing_surface[0]
+    surface = small_grid_and_missing_surface[1]
+    old_fuddy_duddy_crs = rqc.Crs(surface.model, xy_units = 'ft', z_units = 'chain')
+    old_fuddy_duddy_crs.create_xml()
+    surface.model.store_epc()
+    s2 = rqs.Surface(surface.model, uuid = surface.uuid)
+    s2.change_crs(old_fuddy_duddy_crs)
+    name = "test"
+    assert grid.is_aligned
+
+    # Act
+    gcs_optimised = rqgs.find_faces_to_represent_surface_regular_optimised(grid, surface, name)
+
+    # Assert – quite harsh as gcs face ordering could legitimately vary
+    assert gcs_optimised is None