Populate face_face_connectivity (#829)

* Added `face_face_connectivity` construction * Updated Testing Method * add asv --------- Co-authored-by: Philip Chmielowiec <[email protected]> Co-authored-by: Philip Chmielowiec <[email protected]>
UXARRAY · Jul 17, 2024 · 49fd637 · 49fd637
1 parent 0cbbc92
commit 49fd637
Show file tree

Hide file tree

Showing 5 changed files with 65 additions and 4 deletions.
diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py
@@ -92,12 +92,16 @@ class ConnectivityConstruction:
     def setup(self, resolution):
         self.uxds = ux.open_dataset(file_path_dict[resolution][0], file_path_dict[resolution][1])
 
+
     def teardown(self, resolution):
         del self.uxds
 
     def time_n_nodes_per_face(self, resolution):
         self.uxds.uxgrid.n_nodes_per_face
 
+    def time_face_face_connectivity(self, resolution):
+        ux.grid.connectivity._populate_face_face_connectivity(self.uxds.uxgrid)
+
 
 class MatplotlibConversion:
     param_names = ['resolution', 'periodic_elements']

diff --git a/ci/environment.yml b/ci/environment.yml
@@ -32,3 +32,4 @@ dependencies:
   - pip:
     - git+https://github.com/airspeed-velocity/asv
     - pyfma
+    - git+https://github.com/airspeed-velocity/asv
diff --git a/test/test_grid.py b/test/test_grid.py
@@ -8,7 +8,8 @@
 
 import uxarray as ux
 
-from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity, _build_edge_node_connectivity
+from uxarray.grid.connectivity import _populate_face_edge_connectivity, _build_edge_face_connectivity, \
+    _build_edge_node_connectivity, _build_face_face_connectivity, _populate_face_face_connectivity
 
 from uxarray.grid.coordinates import _populate_node_latlon, _lonlat_rad_to_xyz
 
@@ -874,6 +875,23 @@ def test_edge_face_connectivity_sample(self):
         # no invalid entries should occur
         assert n_invalid == 0
 
+    def test_face_face_connectivity_construction(self):
+        """Tests the construction of face-face connectivity."""
+
+        # Open MPAS grid and read in face_face_connectivity
+        grid = ux.open_grid(gridfile_mpas)
+        face_face_conn_old = grid.face_face_connectivity.values
+
+        # Construct new face_face_connectivity using UXarray
+        face_face_conn_new = _build_face_face_connectivity(grid)
+
+        # Sort the arrays before comparison
+        face_face_conn_old_sorted = np.sort(face_face_conn_old, axis=None)
+        face_face_conn_new_sorted = np.sort(face_face_conn_new, axis=None)
+
+        # Assert the new and old face_face_connectivity contains the same faces
+        nt.assert_array_equal(face_face_conn_new_sorted, face_face_conn_old_sorted)
+
 
 class TestClassMethods(TestCase):
     gridfile_ugrid = current_path / "meshfiles" / "ugrid" / "geoflow-small" / "grid.nc"

diff --git a/uxarray/grid/connectivity.py b/uxarray/grid/connectivity.py
@@ -434,3 +434,42 @@ def get_face_node_partitions(n_nodes_per_face):
     change_ind = np.concatenate((np.array([0]), change_ind))
 
     return change_ind, n_nodes_per_face_sorted_ind, element_sizes, size_counts
+
+
+def _populate_face_face_connectivity(grid):
+    """Constructs the UGRID connectivity variable (``face_face_connectivity``)
+    and stores it within the internal (``Grid._ds``) and through the attribute
+    (``Grid.face_face_connectivity``)."""
+    face_face = _build_face_face_connectivity(grid)
+
+    grid._ds["face_face_connectivity"] = xr.DataArray(
+        data=face_face,
+        dims=ugrid.FACE_FACE_CONNECTIVITY_DIMS,
+        attrs=ugrid.FACE_FACE_CONNECTIVITY_ATTRS,
+    )
+
+
+def _build_face_face_connectivity(grid):
+    """Returns face-face connectivity."""
+
+    # Dictionary to store each faces adjacent faces
+    face_neighbors = {i: [] for i in range(grid.n_face)}
+
+    # Loop through each edge_face and add to the dictionary every face that shares an edge
+    for edge_face in grid.edge_face_connectivity.values:
+        face1, face2 = edge_face
+        if face1 != INT_FILL_VALUE and face2 != INT_FILL_VALUE:
+            # Append to each face's dictionary index the opposite face index
+            face_neighbors[face1].append(face2)
+            face_neighbors[face2].append(face1)
+
+    # Convert to an array and pad it with fill values
+    face_face_conn = list(face_neighbors.values())
+    face_face_connectivity = [
+        np.pad(
+            arr, (0, grid.n_max_face_edges - len(arr)), constant_values=INT_FILL_VALUE
+        )
+        for arr in face_face_conn
+    ]
+
+    return face_face_connectivity
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
@@ -38,6 +38,7 @@
     _populate_n_nodes_per_face,
     _populate_node_face_connectivity,
     _populate_edge_face_connectivity,
+    _populate_face_face_connectivity,
 )
 
 from uxarray.grid.geometry import (
@@ -803,9 +804,7 @@ def face_face_connectivity(self) -> xr.DataArray:
         Dimensions ``(n_face, n_max_face_faces)``
         """
         if "face_face_connectivity" not in self._ds:
-            raise NotImplementedError(
-                "Construction of `face_face_connectivity` not yet supported."
-            )
+            _populate_face_face_connectivity(self)
 
         return self._ds["face_face_connectivity"]