From 49fd637bb96b193e9e1d011396f56d4f35a73267 Mon Sep 17 00:00:00 2001 From: Aaron Zedwick <95507181+aaronzedwick@users.noreply.github.com> Date: Wed, 17 Jul 2024 16:54:27 -0500 Subject: [PATCH] Populate `face_face_connectivity` (#829) * Added `face_face_connectivity` construction * Updated Testing Method * add asv --------- Co-authored-by: Philip Chmielowiec <67855069+philipc2@users.noreply.github.com> Co-authored-by: Philip Chmielowiec --- benchmarks/mpas_ocean.py | 4 ++++ ci/environment.yml | 1 + test/test_grid.py | 20 +++++++++++++++++- uxarray/grid/connectivity.py | 39 ++++++++++++++++++++++++++++++++++++ uxarray/grid/grid.py | 5 ++--- 5 files changed, 65 insertions(+), 4 deletions(-) diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py index 71972c5eb..cd1840d82 100644 --- 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 index 9501c397e..88350b007 100644 --- 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 index 01d37a0ae..10faba696 100644 --- 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 index 81a6119b1..41e4201ee 100644 --- 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 index c9b4577a1..ae480b132 100644 --- 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"]