diff --git a/benchmarks/mpas_ocean.py b/benchmarks/mpas_ocean.py
index d53c7165b..468d5779d 100644
--- a/benchmarks/mpas_ocean.py
+++ b/benchmarks/mpas_ocean.py
@@ -185,8 +185,24 @@ def time_const_lat(self, resolution, lat_step):
             self.uxgrid.cross_section.constant_latitude(lat)
 
 
-class PointInPolygon(GridBenchmark):
+class PointInPolygon:
+    param_names = ['resolution']
+    params = ['480km', '120km']
+
+    def setup(self, resolution):
+        self.uxgrid = ux.open_grid(file_path_dict[resolution][0])
+        self.uxgrid.normalize_cartesian_coordinates()
+
+        _ = self.uxgrid.face_edge_connectivity
+
+        point = np.array([0.0, 0.0, 1.0])
+        res = self.uxgrid.get_faces_containing_point(point)
+
+    def teardown(self, resolution):
+        del self.uxgrid
+
     def time_whole_grid(self, resolution):
-        point_xyz = np.array([self.uxgrid.face_x[0].values, self.uxgrid.face_y[0].values, self.uxgrid.face_z[0].values])
+        for i in range(len(self.uxgrid.face_x.values)):
+            point_xyz = np.array([self.uxgrid.face_x[i].values, self.uxgrid.face_y[i].values, self.uxgrid.face_z[i].values])
 
-        self.uxgrid.get_faces_containing_point(point_xyz=point_xyz)
+            self.uxgrid.get_faces_containing_point(point_xyz=point_xyz)
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 487cd44e2..548a2c37d 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -1555,6 +1555,23 @@ def max_face_radius(self):
             self._ds["max_face_radius"] = get_max_face_radius(self)
         return self._ds["max_face_radius"]
 
+    @property
+    def face_edge_nodes_xyz(self):
+        if "face_edge_nodes_xyz" not in self._ds:
+            face_edge = _get_cartesian_face_edge_nodes(
+                self.face_node_connectivity.values,
+                self.n_face,
+                self.n_max_face_nodes,
+                self.node_x.values,
+                self.node_y.values,
+                self.node_z.values,
+            )
+            self._ds["face_edge_nodes_xyz"] = (
+                [self.n_face, self.n_max_face_edges, 2, 3],
+                face_edge,
+            )
+        return self._ds["face_edge_nodes_xyz"]
+
     def chunk(self, n_node="auto", n_edge="auto", n_face="auto"):
         """Converts all arrays to dask arrays with given chunks across grid
         dimensions in-place.
@@ -2432,6 +2449,7 @@ def get_faces_containing_point(self, point_xyz):
 
         # Get the maximum face radius of the grid
         max_face_radius = self.max_face_radius
+        _ = self.face_edge_nodes_xyz
 
         subset = self.subset.bounding_circle(
             center_coord=[*point_xyz],
@@ -2440,19 +2458,11 @@ def get_faces_containing_point(self, point_xyz):
             inverse_indices=True,
         )
 
-        # Get the face's edges for the whole subset
-        face_edge_cartesian = _get_cartesian_face_edge_nodes(
-            subset.face_node_connectivity.values,
-            subset.n_face,
-            subset.n_max_face_nodes,
-            subset.node_x.values,
-            subset.node_y.values,
-            subset.node_z.values,
-        )
+        face_edge_nodes_xyz = subset.face_edge_nodes_xyz.values
 
         inverse_indices = subset.inverse_indices.face.values
 
         # Check if any of the faces in the subset contain the point
-        index = _find_faces(face_edge_cartesian, point_xyz, inverse_indices)
+        index = _find_faces(face_edge_nodes_xyz, point_xyz, inverse_indices)
 
         return index