diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 14709b1a9..eb0371667 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -4,7 +4,6 @@
 
 from html import escape
 
-from numba import njit
 from xarray.core.options import OPTIONS
 
 from typing import (
@@ -14,7 +13,7 @@
     Set,
 )
 
-from uxarray.grid.utils import _get_cartesian_face_edge_nodes
+from uxarray.grid.utils import _get_cartesian_face_edge_nodes, _find_faces
 
 # reader and writer imports
 from uxarray.io._exodus import _read_exodus, _encode_exodus
@@ -71,7 +70,6 @@
     _construct_boundary_edge_indices,
     compute_temp_latlon_array,
     calculate_max_face_radius,
-    point_in_face,
 )
 
 from uxarray.grid.neighbors import (
@@ -2442,6 +2440,7 @@ def get_faces_containing_point(self, point_xyz):
             subset.node_y.values,
             subset.node_z.values,
         )
+
         inverse_indices = subset.inverse_indices.face.values
 
         # Check if any of the faces in the subset contain the point
@@ -2467,25 +2466,3 @@ def get_max_face_radius(self):
         )
 
         return max_distance
-
-
-@njit(cache=True)
-def _find_faces(face_edge_cartesian, point_xyz, inverse_indices):
-    """Finds the faces that contain a given point, inside a subset "face_edge_cartesian"""
-
-    index = []
-
-    # Loop through the whole subset
-    for i, face in enumerate(face_edge_cartesian):
-        # Check if the point is inside the face
-        contains_point = point_in_face(
-            face,
-            point_xyz,
-            inclusive=True,
-        )
-
-        # If the point is, add it to the list of faces
-        if contains_point:
-            index.append(inverse_indices[i])
-
-    return index
diff --git a/uxarray/grid/utils.py b/uxarray/grid/utils.py
index 228fd96fb..11739a01c 100644
--- a/uxarray/grid/utils.py
+++ b/uxarray/grid/utils.py
@@ -5,6 +5,8 @@
 
 from numba import njit
 
+from uxarray.grid.geometry import point_in_face
+
 
 @njit(cache=True)
 def _small_angle_of_2_vectors(u, v):
@@ -458,3 +460,25 @@ def _get_lonlat_rad_face_edge_nodes(
     face_edges_lonlat_rad[valid_mask, 1] = node_lat_rad[valid_edges]
 
     return face_edges_lonlat_rad.reshape(n_face, n_max_face_edges, 2, 2)
+
+
+@njit(cache=True)
+def _find_faces(face_edge_cartesian, point_xyz, inverse_indices):
+    """Finds the faces that contain a given point, inside a subset "face_edge_cartesian"""
+
+    index = []
+
+    # Loop through the whole subset
+    for i, face in enumerate(face_edge_cartesian):
+        # Check if the point is inside the face
+        contains_point = point_in_face(
+            face,
+            point_xyz,
+            inclusive=True,
+        )
+
+        # If the point is, add it to the list of faces
+        if contains_point:
+            index.append(inverse_indices[i])
+
+    return index