Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Normalization of Parsed Cartesian Coordinates #878

Merged
merged 17 commits into from
Sep 18, 2024
Merged
34 changes: 34 additions & 0 deletions uxarray/grid/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@
_set_desired_longitude_range,
_populate_node_latlon,
_populate_node_xyz,
_normalize_xyz,
)
from uxarray.grid.connectivity import (
_populate_edge_node_connectivity,
Expand Down Expand Up @@ -70,6 +71,7 @@
_check_connectivity,
_check_duplicate_nodes,
_check_area,
_check_normalization,
)

from xarray.core.utils import UncachedAccessor
Expand Down Expand Up @@ -1158,6 +1160,38 @@ def compute_face_areas(

return self._face_areas, self._face_jacobian

def normalize_cartesian_coordinates(self):
"""Normalizes Cartesian coordinates."""

if _check_normalization(self):
# check if coordinates are already normalized
return

if "node_x" in self._ds:
# normalize node coordinates
node_x, node_y, node_z = _normalize_xyz(
self.node_x.values, self.node_y.values, self.node_z.values
)
self.node_x.data = node_x
self.node_y.data = node_y
self.node_z.data = node_z
if "edge_x" in self._ds:
# normalize edge coordinates
edge_x, edge_y, edge_z = _normalize_xyz(
self.edge_x.values, self.edge_y.values, self.edge_z.values
)
self.edge_x.data = edge_x
self.edge_y.data = edge_y
self.edge_z.data = edge_z
if "face_x" in self._ds:
# normalize face coordinates
face_x, face_y, face_z = _normalize_xyz(
self.face_x.values, self.face_y.values, self.face_z.values
)
self.face_x.data = face_x
self.face_y.data = face_y
self.face_z.data = face_z

def to_xarray(self, grid_format: Optional[str] = "ugrid"):
"""Returns a xarray Dataset representation in a specific grid format
from the Grid object.
Expand Down
32 changes: 24 additions & 8 deletions uxarray/grid/validation.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@


# validation helper functions
def _check_connectivity(self):
def _check_connectivity(grid):
"""Check if all nodes are referenced by at least one element.

If not, the mesh may have hanging nodes and may not a valid UGRID
Expand All @@ -14,28 +14,28 @@ def _check_connectivity(self):

# Check if all nodes are referenced by at least one element
# get unique nodes in connectivity
nodes_in_conn = np.unique(self.face_node_connectivity.values.flatten())
nodes_in_conn = np.unique(grid.face_node_connectivity.values.flatten())
# remove negative indices/fill values from the list
nodes_in_conn = nodes_in_conn[nodes_in_conn >= 0]

# check if the size of unique nodes in connectivity is equal to the number of nodes
if nodes_in_conn.size == self.n_node:
if nodes_in_conn.size == grid.n_node:
print("-All nodes are referenced by at least one element.")
return True
else:
warn(
"Some nodes may not be referenced by any element. {0} and {1}".format(
nodes_in_conn.size, self.n_node
nodes_in_conn.size, grid.n_node
),
RuntimeWarning,
)
return False


def _check_duplicate_nodes(self):
def _check_duplicate_nodes(grid):
"""Check if there are duplicate nodes in the mesh."""

coords1 = np.column_stack((np.vstack(self.node_lon), np.vstack(self.node_lat)))
coords1 = np.column_stack((np.vstack(grid.node_lon), np.vstack(grid.node_lat)))
unique_nodes, indices = np.unique(coords1, axis=0, return_index=True)
duplicate_indices = np.setdiff1d(np.arange(len(coords1)), indices)

Expand All @@ -52,9 +52,9 @@ def _check_duplicate_nodes(self):
return True


def _check_area(self):
def _check_area(grid):
"""Check if each face area is greater than our constant ERROR_TOLERANCE."""
areas = self.face_areas
areas = grid.face_areas
# Check if area of any face is close to zero
if np.any(np.isclose(areas, 0, atol=ERROR_TOLERANCE)):
warn(
Expand All @@ -65,3 +65,19 @@ def _check_area(self):
else:
print("-No face area is close to zero.")
return True


def _check_normalization(grid):
"""Checks whether all the cartesiain coordinates are normalized."""

if "node_x" in grid._ds:
if not (grid.node_x**2 + grid.node_y**2 + grid.node_z**2).all():
return False
if "edge_x" in grid._ds:
if not (grid.edge_x**2 + grid.edge_y**2 + grid.edge_z**2).all():
return False
if "face_x" in grid._ds:
if not (grid.face_x**2 + grid.face_y**2 + grid.face_z**2).all():
return False

return True
philipc2 marked this conversation as resolved.
Show resolved Hide resolved
Loading