diff --git a/docs/changelog.rst b/docs/changelog.rst
index c364491dd..62c793d71 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -6,6 +6,27 @@ All notable changes to this project will be documented in this file.
 The format is based on `Keep a Changelog`_, and this project adheres to
 `Semantic Versioning`_.
 
+Unreleased
+----------
+
+Fixed
+~~~~~
+
+- Release 0.12.0 changed the return type of the face node connectivity of
+  :attr:`xugrid.Ugrid2d.voronoi_topology` from a `scipy.sparse.coo_matrix` to
+  an ordinary `np.array` of integers (and similarly for internal voronoi
+  tesselations); this dense array had fill (hard-coded) values of -1,
+  potentially differing from the grid's fill value. This lead to a number of
+  errors for methods relying on voronoi tesselations (such as contour plots)
+  if the fill value of the grid was not -1. Internally, a ``FILL_VALUE = -1``
+  is now used everywhere in connectivity arrays, and fill values are no longer
+  passed for internal methods; a value of -1 is always assumed. When converting
+  the grid (back) to a dataset with :meth:`xugrid.Ugrid1d.to_dataset` or
+  :meth:`xugrid.Ugrid2d.to_dataset`, the fill value is set back to its original
+  value; the fill value is also set when calling
+  :meth:`xugrid.UgridDataArrayAccessor.to_netcdf` or
+  :meth:`xugrid.UgridDatasetAccessor.to_netcdf`.
+
 [0.12.0] 2024-09-03
 -------------------
 
diff --git a/examples-dev/voronoi.py b/examples-dev/voronoi.py
index 758ea5a70..cd5779e9e 100644
--- a/examples-dev/voronoi.py
+++ b/examples-dev/voronoi.py
@@ -19,11 +19,14 @@
 modules, should you not want to rely on more complex dependencies such as
 ``xugrid`` and ``xarray``.
 """
+# %%
 import matplotlib.pyplot as plt
 import matplotlib.tri as mtri
 import numpy as np
 from matplotlib.collections import LineCollection, PolyCollection
 
+from xugrid.constants import FILL_VALUE  # equals -1
+
 # %%
 # From xugrid, we need only import the ``connectivity`` and ``voronoi``
 # modules. The functions in these modules depend only on ``numpy`` and
@@ -63,12 +66,12 @@ def generate_disk(partitions: int, depth: int):
     return np.column_stack((x, y)), triang.triangles
 
 
-def edge_plot(vertices, edge_nodes, ax, fill_value=-1, **kwargs):
+def edge_plot(vertices, edge_nodes, ax, **kwargs):
     n_edge = len(edge_nodes)
     edge_coords = np.empty((n_edge, 2, 2), dtype=float)
     node_0 = edge_nodes[:, 0]
     node_1 = edge_nodes[:, 1]
-    valid = (node_0 != fill_value) & (node_1 != fill_value)
+    valid = (node_0 != FILL_VALUE) & (node_1 != FILL_VALUE)
     node_0 = node_0[valid]
     node_1 = node_1[valid]
     edge_coords[:, 0, 0] = vertices[node_0, 0]
@@ -81,10 +84,10 @@ def edge_plot(vertices, edge_nodes, ax, fill_value=-1, **kwargs):
     return primitive
 
 
-def face_plot(vertices, face_nodes, ax, fill_value=-1, **kwargs):
+def face_plot(vertices, face_nodes, ax, **kwargs):
     vertices = vertices[face_nodes]
     # Replace fill value; PolyCollection ignores NaN.
-    vertices[face_nodes == fill_value] = np.nan
+    vertices[face_nodes == FILL_VALUE] = np.nan
     collection = PolyCollection(vertices, **kwargs)
     primitive = ax.add_collection(collection)
     ax.autoscale()
@@ -107,12 +110,12 @@ def comparison_plot(
         sharex=True,
     )
 
-    edges0, _ = connectivity.edge_connectivity(faces0, -1)
+    edges0, _ = connectivity.edge_connectivity(faces0)
     edge_plot(vertices0, edges0, ax0, colors="black")
     ax0.scatter(*centroids0.T, color="red")
     ax0.scatter(*vertices0.T, color="black")
 
-    edges1, _ = connectivity.edge_connectivity(faces1, -1)
+    edges1, _ = connectivity.edge_connectivity(faces1)
     edge_plot(vertices0, edges0, ax1, colors="black")
     edge_plot(vertices1, edges1, ax1, colors="red")
 
@@ -126,13 +129,11 @@ def comparison_plot(
 # %%
 # Let's start by generating a simple unstructured mesh and use only its
 # centroids to generate a voronoi tesselation.
-#
-# Note: ``-1`` functions as the fill value in this example.
 
 vertices, faces = generate_disk(5, 2)
 centroids = vertices[faces].mean(axis=1)
 
-node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
@@ -158,17 +159,14 @@ def comparison_plot(
 # The ``voronoi_topology`` is capable of preserving the exterior exactly, but
 # this requires more topological information.
 
-edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
-    faces, -1
-)
-edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(faces)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
 )
@@ -190,18 +188,15 @@ def comparison_plot(
 faces = connectivity.renumber(new)
 centroids = vertices[faces].mean(axis=1)
 
-node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
-edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
-    faces, -1
-)
-edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(faces)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
 )
@@ -217,18 +212,15 @@ def comparison_plot(
 # of the original mesh altogether. We still add an orthogonal projection of
 # every centroid to exterior edges.
 
-node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
-edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
-    faces, -1
-)
-edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(faces)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=False,
 )
@@ -240,18 +232,15 @@ def comparison_plot(
 # exactly. Alternatively, we can choose to skip the exterior vertex if it
 # creates a concave face:
 
-node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
-edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
-    faces, -1
-)
-edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(faces)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
 voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
     skip_concave=True,
@@ -267,7 +256,7 @@ def comparison_plot(
     vertices,
     centroids,
 )
-edges0, _ = connectivity.edge_connectivity(faces0, -1)
+edges0, _ = connectivity.edge_connectivity(faces0)
 
 nodes1, faces1, face_index1, _ = voronoi.voronoi_topology(
     node_face_connectivity,
@@ -275,11 +264,10 @@ def comparison_plot(
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=False,
 )
-edges1, _ = connectivity.edge_connectivity(faces1, -1)
+edges1, _ = connectivity.edge_connectivity(faces1)
 
 nodes2, faces2, _, _ = voronoi.voronoi_topology(
     node_face_connectivity,
@@ -287,11 +275,10 @@ def comparison_plot(
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
 )
-edges2, _ = connectivity.edge_connectivity(faces2, -1)
+edges2, _ = connectivity.edge_connectivity(faces2)
 
 nodes3, faces3, face_index3, node_map3 = voronoi.voronoi_topology(
     node_face_connectivity,
@@ -299,12 +286,11 @@ def comparison_plot(
     centroids,
     edge_face_connectivity=edge_face_connectivity,
     edge_node_connectivity=edge_node_connectivity,
-    fill_value=-1,
     add_exterior=True,
     add_vertices=True,
     skip_concave=True,
 )
-edges3, _ = connectivity.edge_connectivity(faces3, -1)
+edges3, _ = connectivity.edge_connectivity(faces3)
 
 fig, axes = plt.subplots(
     nrows=1,
@@ -342,10 +328,10 @@ def comparison_plot(
 # fourth option, since it includes some vertices of the original mesh, which
 # are connected to multiple faces.
 
-triangles0, face_triangles0 = connectivity.triangulate(faces0, -1)
+triangles0, face_triangles0 = connectivity.triangulate(faces0)
 triangulation0 = mtri.Triangulation(nodes0[:, 0], nodes0[:, 1], triangles0)
 
-triangles1, face_triangles1 = connectivity.triangulate(faces1, -1)
+triangles1, face_triangles1 = connectivity.triangulate(faces1)
 triangulation1 = mtri.Triangulation(nodes1[:, 0], nodes1[:, 1], triangles1)
 
 
@@ -422,3 +408,5 @@ def comparison_plot(
 
 ax0.set_xlim(-1.5, 1.5)
 ax0.set_ylim(-1.5, 1.5)
+
+# %%
diff --git a/tests/test_connectivity.py b/tests/test_connectivity.py
index 1c980c8cc..9be9bf636 100644
--- a/tests/test_connectivity.py
+++ b/tests/test_connectivity.py
@@ -2,12 +2,12 @@
 import pytest
 from scipy import sparse
 
+from xugrid.constants import FILL_VALUE
 from xugrid.ugrid import connectivity
 
 
 @pytest.fixture(scope="function")
 def triangle_mesh():
-    fill_value = -1
     # Two triangles
     faces = np.array(
         [
@@ -15,20 +15,19 @@ def triangle_mesh():
             [1, 3, 2],
         ]
     )
-    return faces, fill_value
+    return faces
 
 
 @pytest.fixture(scope="function")
 def mixed_mesh():
-    fill_value = -1
     # Triangle, quadrangle
     faces = np.array(
         [
-            [0, 1, 2, fill_value],
+            [0, 1, 2, FILL_VALUE],
             [1, 3, 4, 2],
         ]
     )
-    return faces, fill_value
+    return faces
 
 
 def test_argsort_rows():
@@ -88,39 +87,39 @@ def test_neighbors():
 
 
 def test_to_ij(triangle_mesh, mixed_mesh):
-    faces, fill_value = triangle_mesh
-    actual_i, actual_j = connectivity._to_ij(faces, fill_value, invert=False)
+    faces = triangle_mesh
+    actual_i, actual_j = connectivity._to_ij(faces, invert=False)
     expected_i = [0, 0, 0, 1, 1, 1]
     expected_j = [0, 1, 2, 1, 3, 2]
     assert np.array_equal(actual_i, expected_i)
     assert np.array_equal(actual_j, expected_j)
 
     # Inverted
-    actual_i, actual_j = connectivity._to_ij(faces, fill_value, invert=True)
+    actual_i, actual_j = connectivity._to_ij(faces, invert=True)
     assert np.array_equal(actual_i, expected_j)
     assert np.array_equal(actual_j, expected_i)
 
-    faces, fill_value = mixed_mesh
-    actual_i, actual_j = connectivity._to_ij(faces, fill_value, invert=False)
+    faces = mixed_mesh
+    actual_i, actual_j = connectivity._to_ij(faces, invert=False)
     expected_i = [0, 0, 0, 1, 1, 1, 1]
     expected_j = [0, 1, 2, 1, 3, 4, 2]
     assert np.array_equal(actual_i, expected_i)
     assert np.array_equal(actual_j, expected_j)
 
     # Inverted
-    actual_i, actual_j = connectivity._to_ij(faces, fill_value, invert=True)
+    actual_i, actual_j = connectivity._to_ij(faces, invert=True)
     assert np.array_equal(actual_i, expected_j)
     assert np.array_equal(actual_j, expected_i)
 
 
 def test_to_sparse(mixed_mesh):
-    faces, fill_value = mixed_mesh
-    csr = connectivity._to_sparse(faces, fill_value, invert=False, sort_indices=True)
+    faces = mixed_mesh
+    csr = connectivity._to_sparse(faces, invert=False, sort_indices=True)
     expected_j = np.array([0, 1, 2, 1, 2, 3, 4])
     assert np.array_equal(csr.indices, expected_j)
     assert csr.has_sorted_indices
 
-    csr = connectivity._to_sparse(faces, fill_value, invert=False, sort_indices=False)
+    csr = connectivity._to_sparse(faces, invert=False, sort_indices=False)
     expected_j = np.array([0, 1, 2, 1, 3, 4, 2])
     assert np.array_equal(csr.indices, expected_j)
     assert not csr.has_sorted_indices
@@ -142,23 +141,23 @@ def test_ragged_index():
 
 
 def test_sparse_dense_conversion_roundtrip(triangle_mesh, mixed_mesh):
-    faces, fill_value = triangle_mesh
-    sparse = connectivity.to_sparse(faces, fill_value)
-    back = connectivity.to_dense(sparse, fill_value)
+    faces = triangle_mesh
+    sparse = connectivity.to_sparse(faces)
+    back = connectivity.to_dense(sparse)
     # Note: roundtrip does not preserve CW/CCW orientation, since orientation
     # does not apply to node_face_connectivity, but the sorted rows should
     # contain the same elements.
     assert np.array_equal(faces.sort(axis=1), back.sort(axis=1))
 
-    faces, fill_value = mixed_mesh
-    sparse = connectivity.to_sparse(faces, fill_value)
-    back = connectivity.to_dense(sparse, fill_value)
+    faces = mixed_mesh
+    sparse = connectivity.to_sparse(faces)
+    back = connectivity.to_dense(sparse)
     assert np.array_equal(faces.sort(axis=1), back.sort(axis=1))
 
 
 def test_invert_dense(triangle_mesh, mixed_mesh):
-    faces, fill_value = triangle_mesh
-    actual = connectivity.invert_dense(faces, fill_value)
+    faces = triangle_mesh
+    actual = connectivity.invert_dense(faces)
     expected = np.array(
         [
             [0, -1],  # 0
@@ -169,8 +168,8 @@ def test_invert_dense(triangle_mesh, mixed_mesh):
     )
     assert np.array_equal(actual, expected)
 
-    faces, fill_value = mixed_mesh
-    actual = connectivity.invert_dense(faces, fill_value)
+    faces = mixed_mesh
+    actual = connectivity.invert_dense(faces)
     expected = np.array(
         [
             [0, -1],  # 0
@@ -184,10 +183,10 @@ def test_invert_dense(triangle_mesh, mixed_mesh):
 
 
 def test_invert_sparse(triangle_mesh, mixed_mesh):
-    faces, fill_value = triangle_mesh
-    sparse = connectivity.to_sparse(faces, fill_value)
+    faces = triangle_mesh
+    sparse = connectivity.to_sparse(faces)
     inverted = connectivity.invert_sparse(sparse)
-    actual = connectivity.to_dense(inverted, fill_value)
+    actual = connectivity.to_dense(inverted)
     expected = np.array(
         [
             [0, -1],  # 0
@@ -198,10 +197,10 @@ def test_invert_sparse(triangle_mesh, mixed_mesh):
     )
     assert np.array_equal(actual, expected)
 
-    faces, fill_value = mixed_mesh
-    sparse = connectivity.to_sparse(faces, fill_value)
+    faces = mixed_mesh
+    sparse = connectivity.to_sparse(faces)
     inverted = connectivity.invert_sparse(sparse)
-    actual = connectivity.to_dense(inverted, fill_value)
+    actual = connectivity.to_dense(inverted)
     expected = np.array(
         [
             [0, -1],  # 0
@@ -215,30 +214,30 @@ def test_invert_sparse(triangle_mesh, mixed_mesh):
 
 
 def test_to_dense(triangle_mesh):
-    faces, fill_value = triangle_mesh
-    sparse = connectivity.to_sparse(faces, fill_value)
-    actual = connectivity.to_dense(sparse, fill_value)
+    faces = triangle_mesh
+    sparse = connectivity.to_sparse(faces)
+    actual = connectivity.to_dense(sparse)
     assert np.array_equal(actual, np.sort(faces, axis=1))
 
     with pytest.raises(ValueError, match="n_columns 2 is too small"):
-        connectivity.to_dense(sparse, fill_value, n_columns=2)
+        connectivity.to_dense(sparse, n_columns=2)
 
     # now pad
-    actual = connectivity.to_dense(sparse, fill_value, n_columns=4)
+    actual = connectivity.to_dense(sparse, n_columns=4)
     expected = np.array(
         [
-            [0, 1, 2, fill_value],
-            [1, 2, 3, fill_value],
+            [0, 1, 2, FILL_VALUE],
+            [1, 2, 3, FILL_VALUE],
         ]
     )
     assert np.array_equal(actual, expected)
 
     # and twice
-    actual = connectivity.to_dense(sparse, fill_value, n_columns=5)
+    actual = connectivity.to_dense(sparse, n_columns=5)
     expected = np.array(
         [
-            [0, 1, 2, fill_value, fill_value],
-            [1, 2, 3, fill_value, fill_value],
+            [0, 1, 2, FILL_VALUE, FILL_VALUE],
+            [1, 2, 3, FILL_VALUE, FILL_VALUE],
         ]
     )
     assert np.array_equal(actual, expected)
@@ -283,12 +282,12 @@ def test_renumber():
 def test_renumber_with_fill_value():
     a = np.array(
         [
-            [0, 1, -1],
+            [0, 1, FILL_VALUE],
             [10, 11, 12],
-            [30, -1, 32],
+            [30, FILL_VALUE, 32],
         ]
     )
-    actual = connectivity.renumber(a, fill_value=-1)
+    actual = connectivity.renumber(a)
     expected = np.array(
         [
             [0, 1, -1],
@@ -305,7 +304,7 @@ def test_renumber_with_fill_value():
             [30, -1, 2],
         ]
     )
-    actual = connectivity.renumber(a, fill_value=-1)
+    actual = connectivity.renumber(a)
     expected = np.array(
         [
             [0, -1, 1],
@@ -317,8 +316,8 @@ def test_renumber_with_fill_value():
 
 
 def test_close_polygons(mixed_mesh):
-    faces, fill_value = mixed_mesh
-    closed, isfill = connectivity.close_polygons(faces, fill_value)
+    faces = mixed_mesh
+    closed, isfill = connectivity.close_polygons(faces)
     expected = np.array(
         [
             [0, 1, 2, 0, 0],
@@ -333,11 +332,11 @@ def test_close_polygons(mixed_mesh):
 
 
 def test_reverse_orientation(mixed_mesh):
-    faces, fill_value = mixed_mesh
-    reverse = connectivity.reverse_orientation(faces, fill_value)
+    faces = mixed_mesh
+    reverse = connectivity.reverse_orientation(faces)
     expected = np.array(
         [
-            [2, 1, 0, fill_value],
+            [2, 1, 0, FILL_VALUE],
             [2, 4, 3, 1],
         ]
     )
@@ -353,32 +352,31 @@ def test_counterclockwise():
             [0.0, 2.0],
         ]
     )
-    fill_value = -1
 
     # Already ccw, nothing should be changed.
     faces = np.array([[0, 2, 3, -1]])
-    actual = connectivity.counterclockwise(faces, fill_value, nodes)
+    actual = connectivity.counterclockwise(faces, nodes)
     assert np.array_equal(actual, faces)
 
     # Clockwise with a fill value, reverse.
     faces_cw = np.array([[3, 2, 0, -1]])
-    actual = connectivity.counterclockwise(faces_cw, fill_value, nodes)
+    actual = connectivity.counterclockwise(faces_cw, nodes)
     assert np.array_equal(actual, faces)
 
     # Including a hanging node, ccw, nothing changed.
     hanging_ccw = np.array([[0, 1, 2, 3, -1]])
-    actual = connectivity.counterclockwise(hanging_ccw, fill_value, nodes)
+    actual = connectivity.counterclockwise(hanging_ccw, nodes)
     assert np.array_equal(actual, hanging_ccw)
 
     # Including a hanging node, reverse.
     hanging_cw = np.array([[3, 2, 1, 0, -1]])
-    actual = connectivity.counterclockwise(hanging_cw, fill_value, nodes)
+    actual = connectivity.counterclockwise(hanging_cw, nodes)
     assert np.array_equal(actual, hanging_ccw)
 
 
 def test_edge_connectivity(mixed_mesh):
-    faces, fill_value = mixed_mesh
-    edge_nodes, face_edges = connectivity.edge_connectivity(faces, fill_value)
+    faces = mixed_mesh
+    edge_nodes, face_edges = connectivity.edge_connectivity(faces)
     expected_edge_nodes = np.array(
         [
             [0, 1],
@@ -400,10 +398,10 @@ def test_edge_connectivity(mixed_mesh):
 
 
 def test_validate_edge_connectivity(mixed_mesh):
-    faces, fill_value = mixed_mesh
+    faces = mixed_mesh
     edges = np.array([[0, 1]])
     with pytest.raises(ValueError, match="face_node_connectivity defines 6 edges"):
-        connectivity.validate_edge_node_connectivity(faces, fill_value, edges)
+        connectivity.validate_edge_node_connectivity(faces, edges)
 
     edges = np.array(
         [
@@ -418,7 +416,7 @@ def test_validate_edge_connectivity(mixed_mesh):
             [0, 4],  # F
         ]
     )
-    actual = connectivity.validate_edge_node_connectivity(faces, fill_value, edges)
+    actual = connectivity.validate_edge_node_connectivity(faces, edges)
     expected = np.array([True, False, False, True, True, True, True, True, False])
     assert np.array_equal(actual, expected)
 
@@ -505,7 +503,7 @@ def test_face_face_connectivity():
             [1, -1],
         ]
     )
-    face_face = connectivity.face_face_connectivity(edge_faces, fill_value=-1)
+    face_face = connectivity.face_face_connectivity(edge_faces)
     assert isinstance(face_face, sparse.csr_matrix)
     assert np.array_equal(face_face.indices, [1, 0])
     assert np.array_equal(face_face.indptr, [0, 1, 2])
@@ -513,7 +511,7 @@ def test_face_face_connectivity():
 
 
 def test_centroids(mixed_mesh):
-    faces, fill_value = mixed_mesh
+    faces = mixed_mesh
     nodes = np.array(
         [
             [0.0, 0.0],
@@ -523,7 +521,7 @@ def test_centroids(mixed_mesh):
             [2.0, 1.0],
         ]
     )
-    actual = connectivity.centroids(faces, fill_value, nodes[:, 0], nodes[:, 1])
+    actual = connectivity.centroids(faces, nodes[:, 0], nodes[:, 1])
     expected = np.array(
         [
             [2.0 / 3.0, 1.0 / 3.0],
@@ -534,7 +532,7 @@ def test_centroids(mixed_mesh):
 
 
 def test_circumcenters_error(mixed_mesh):
-    faces, fill_value = mixed_mesh
+    faces = mixed_mesh
     nodes = np.array(
         [
             [0.0, 0.0],
@@ -545,11 +543,11 @@ def test_circumcenters_error(mixed_mesh):
         ]
     )
     with pytest.raises(NotImplementedError):
-        connectivity.circumcenters(faces, fill_value, nodes[:, 0], nodes[:, 1])
+        connectivity.circumcenters(faces, nodes[:, 0], nodes[:, 1])
 
 
 def test_circumcenters(triangle_mesh):
-    faces, fill_value = triangle_mesh
+    faces = triangle_mesh
     nodes = np.array(
         [
             [0.0, 1.0],
@@ -558,7 +556,7 @@ def test_circumcenters(triangle_mesh):
             [1.0, -2.0],
         ]
     )
-    actual = connectivity.circumcenters(faces, fill_value, nodes[:, 0], nodes[:, 1])
+    actual = connectivity.circumcenters(faces, nodes[:, 0], nodes[:, 1])
     expected = np.array([[1.0, 0.5], [1.0, -0.75]])
     assert np.allclose(actual, expected)
 
@@ -585,20 +583,20 @@ def test_structured_connectivity():
 def test_area(mixed_mesh):
     node_x = np.array([0.0, 1.0, 1.0, 2.0, 2.0])
     node_y = np.array([0.0, 0.0, 1.0, 0.0, 1.0])
-    actual = connectivity.area(*mixed_mesh, node_x, node_y)
+    actual = connectivity.area(mixed_mesh, node_x, node_y)
     assert np.allclose(actual, [0.5, 1.0])
 
 
 def test_perimeter(mixed_mesh):
     node_x = np.array([0.0, 1.0, 1.0, 2.0, 2.0])
     node_y = np.array([0.0, 0.0, 1.0, 0.0, 1.0])
-    actual = connectivity.perimeter(*mixed_mesh, node_x, node_y)
+    actual = connectivity.perimeter(mixed_mesh, node_x, node_y)
     assert np.allclose(actual, [2.0 + np.sqrt(2), 4.0])
 
 
 def test_triangulate(mixed_mesh):
-    faces, fill_value = mixed_mesh
-    actual_triangles, actual_faces = connectivity.triangulate_dense(faces, fill_value)
+    faces = mixed_mesh
+    actual_triangles, actual_faces = connectivity.triangulate_dense(faces)
     expected_triangles = np.array(
         [
             [0, 1, 2],
@@ -610,7 +608,7 @@ def test_triangulate(mixed_mesh):
     assert np.array_equal(actual_triangles, expected_triangles)
     assert np.array_equal(actual_faces, expected_faces)
 
-    sparse_faces = connectivity.to_sparse(faces, -1, sort_indices=False).tocoo()
+    sparse_faces = connectivity.to_sparse(faces, sort_indices=False).tocoo()
     actual_triangles, actual_faces = connectivity.triangulate_coo(sparse_faces)
     assert np.array_equal(actual_triangles, expected_triangles)
     assert np.array_equal(actual_faces, expected_faces)
diff --git a/tests/test_conversion.py b/tests/test_conversion.py
index 12936314f..abfdf67f5 100644
--- a/tests/test_conversion.py
+++ b/tests/test_conversion.py
@@ -5,6 +5,7 @@
 import xarray as xr
 
 from xugrid import conversion as cv
+from xugrid.constants import FILL_VALUE
 
 
 @pytest.fixture(scope="function")
@@ -32,7 +33,6 @@ def line_gdf():
 def triangle_mesh():
     x = np.array([0.0, 1.0, 1.0, 2.0])
     y = np.array([0.0, 0.0, 1.0, 0.0])
-    fill_value = -1
     # Two triangles
     faces = np.array(
         [
@@ -40,22 +40,21 @@ def triangle_mesh():
             [1, 3, 2],
         ]
     )
-    return x, y, faces, fill_value
+    return x, y, faces
 
 
 @pytest.fixture(scope="function")
 def mixed_mesh():
     x = np.array([0.0, 1.0, 1.0, 2.0, 2.0])
     y = np.array([0.0, 0.0, 1.0, 0.0, 1.0])
-    fill_value = -1
     # Triangle, quadrangle
     faces = np.array(
         [
-            [0, 1, 2, fill_value],
+            [0, 1, 2, FILL_VALUE],
             [1, 3, 4, 2],
         ]
     )
-    return x, y, faces, fill_value
+    return x, y, faces
 
 
 @pytest.fixture(scope="function")
@@ -111,14 +110,13 @@ def test_edges_geos_roundtrip(line):
 # Cannot use fixtures in parametrize:
 # https://github.com/pytest-dev/pytest/issues/349
 def _faces_geos_roundtrip(mesh):
-    x, y, c, fv = mesh
-    actual = cv.faces_to_polygons(x, y, c, fv)
-    x_back, y_back, c_back, fv_back = cv.polygons_to_faces(actual)
-    polygons_back = cv.faces_to_polygons(x_back, y_back, c_back, fv_back)
+    x, y, c = mesh
+    actual = cv.faces_to_polygons(x, y, c)
+    x_back, y_back, c_back = cv.polygons_to_faces(actual)
+    polygons_back = cv.faces_to_polygons(x_back, y_back, c_back)
     assert np.array_equal(x, x_back)
     assert np.array_equal(y, y_back)
     assert np.array_equal(c, c_back)
-    assert fv == fv_back
     assert np.array_equal(actual, polygons_back)
 
 
diff --git a/tests/test_ugrid2d.py b/tests/test_ugrid2d.py
index ee787f93a..944816a88 100644
--- a/tests/test_ugrid2d.py
+++ b/tests/test_ugrid2d.py
@@ -1051,6 +1051,19 @@ def test_tesselate_centroidal_voronoi():
     voronoi = grid.tesselate_centroidal_voronoi()
     assert voronoi.n_face == 7
 
+    faces = FACES.copy()
+    faces[faces == -1] = -999
+    grid = xugrid.Ugrid2d(
+        node_x=VERTICES[:, 0],
+        node_y=VERTICES[:, 1],
+        fill_value=-999,
+        face_node_connectivity=faces,
+    )
+    voronoi = grid.tesselate_centroidal_voronoi(add_exterior=True)
+    vfaces = voronoi.face_node_connectivity
+    fill_nodes = vfaces[vfaces < 0]
+    assert (fill_nodes == -1).all()
+
 
 def test_tesselate_circumcenter_voronoi():
     grid = grid2d()
diff --git a/tests/test_ugrid_dataset.py b/tests/test_ugrid_dataset.py
index 39add6b7b..48e65ced5 100644
--- a/tests/test_ugrid_dataset.py
+++ b/tests/test_ugrid_dataset.py
@@ -1191,6 +1191,20 @@ def test_fm_fillvalue_startindex_isel():
     # xugrid 0.6.0 raises "ValueError: Invalid edge_node_connectivity"
     uds.isel({uds.grid.face_dimension: [1]})
 
+    # Check internal fill value. Should be FILL_VALUE
+    grid = uds.ugrid.grid
+    assert (grid.face_node_connectivity != -999).all()
+    gridds = grid.to_dataset()
+    # Should be set back to the origina fill value.
+    assert (gridds["mesh2d_face_nodes"] != xugrid.constants.FILL_VALUE).all()
+
+    # And similarly for the UgridAccessors.
+    ds = uds.ugrid.to_dataset()
+    assert (ds["mesh2d_face_nodes"] != xugrid.constants.FILL_VALUE).all()
+
+    ds_uda = uds["mesh2d_facevar"].ugrid.to_dataset()
+    assert (ds_uda["mesh2d_face_nodes"] != xugrid.constants.FILL_VALUE).all()
+
 
 def test_fm_facenodeconnectivity_fillvalue():
     """
diff --git a/tests/test_voronoi.py b/tests/test_voronoi.py
index 9c68c92f6..eea8fd5b2 100644
--- a/tests/test_voronoi.py
+++ b/tests/test_voronoi.py
@@ -79,7 +79,6 @@ def setup(self):
                 [3.0, 2.0],  # 11
             ]
         )
-        self.fill_value = -1
         self.face_node_connectivity = np.array(
             [
                 [0, 1, 5, 4],
@@ -91,18 +90,15 @@ def setup(self):
             ]
         )
         self.node_face_connectivity = connectivity.invert_dense_to_sparse(
-            self.face_node_connectivity, self.fill_value
+            self.face_node_connectivity
         )
         (
             self.edge_node_connectivity,
             face_edge_connectivity,
         ) = connectivity.edge_connectivity(
             self.face_node_connectivity,
-            self.fill_value,
-        )
-        self.edge_face_connectivity = connectivity.invert_dense(
-            face_edge_connectivity, self.fill_value
         )
+        self.edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity)
         self.centroids = np.array(
             [
                 [0.5, 0.5],
@@ -160,7 +156,6 @@ def test_interior_centroids(self):
             self.node_face_connectivity,
             self.edge_face_connectivity,
             self.edge_node_connectivity,
-            fill_value=self.fill_value,
         )
         expected_i = np.array([1, 1, 2, 2, 4, 4, 7, 7, 9, 9, 10, 10])
         expected_j = np.array([0, 1, 1, 2, 0, 3, 2, 5, 3, 4, 4, 5])
@@ -178,7 +173,6 @@ def test_exterior_vertices(self):
         ) = voronoi.exterior_vertices(
             self.edge_face_connectivity,
             self.edge_node_connectivity,
-            self.fill_value,
             self.vertices,
             self.centroids,
             add_vertices=False,
@@ -223,7 +217,6 @@ def test_voronoi_topology__add_exterior(self):
             self.centroids,
             self.edge_face_connectivity,
             self.edge_node_connectivity,
-            self.fill_value,
             add_exterior=True,
         )
         expected_vertices = rowsort(
@@ -240,7 +233,6 @@ def test_voronoi_topology__add_exterior(self):
             self.centroids,
             self.edge_face_connectivity,
             self.edge_node_connectivity,
-            self.fill_value,
             add_exterior=True,
             add_vertices=True,
         )
diff --git a/xugrid/constants.py b/xugrid/constants.py
index befb2a387..4df367236 100644
--- a/xugrid/constants.py
+++ b/xugrid/constants.py
@@ -22,6 +22,10 @@
 PolygonArray = np.ndarray
 SparseMatrix = Union[coo_matrix, csr_matrix]
 
+# Internally we always use a fill value of -1. This ensures we can always index
+# with the fill value as well, since any array will have at least size 1.
+FILL_VALUE = -1
+
 
 class Point(NamedTuple):
     x: float
diff --git a/xugrid/conversion.py b/xugrid/conversion.py
index 465fcc46c..9b1e76ec5 100644
--- a/xugrid/conversion.py
+++ b/xugrid/conversion.py
@@ -11,6 +11,7 @@
 import xarray as xr
 
 from xugrid.constants import (
+    FILL_VALUE,
     FloatArray,
     IntArray,
     IntDType,
@@ -48,9 +49,9 @@ def edges_to_linestrings(
 
 
 def faces_to_polygons(
-    x: FloatArray, y: FloatArray, face_node_connectivity: IntArray, fill_value: int
+    x: FloatArray, y: FloatArray, face_node_connectivity: IntArray
 ) -> PolygonArray:
-    is_data = face_node_connectivity != fill_value
+    is_data = face_node_connectivity != FILL_VALUE
     m_per_row = is_data.sum(axis=1)
     i = np.repeat(np.arange(len(face_node_connectivity)), m_per_row)
     c = face_node_connectivity.ravel()[is_data.ravel()]
@@ -98,7 +99,6 @@ def polygons_to_faces(
     n = len(polygons)
     m_per_row = np.bincount(indices)
     m = m_per_row.max()
-    fill_value = -1
     # Allocate 2D array and create a flat view of the dense connectivity
     conn = np.empty((n, m), dtype=IntDType)
     flat_conn = conn.ravel()
@@ -108,10 +108,10 @@ def polygons_to_faces(
         valid = slice(None)  # a[:] equals a[slice(None)]
     else:
         valid = ragged_index(n, m, m_per_row).ravel()
-        flat_conn[~valid] = fill_value
+        flat_conn[~valid] = FILL_VALUE
     flat_conn[valid] = inverse
     x, y = contiguous_xy(unique)
-    return x, y, conn, fill_value
+    return x, y, conn
 
 
 def _scalar_spacing(coords, spacing):
diff --git a/xugrid/plot/plot.py b/xugrid/plot/plot.py
index cc6a425ba..7f573b2b7 100644
--- a/xugrid/plot/plot.py
+++ b/xugrid/plot/plot.py
@@ -511,7 +511,7 @@ def pcolormesh(grid, da, ax, **kwargs):
         raise ValueError("pcolormesh only supports data on faces")
 
     nodes = grid.node_coordinates
-    faces, _ = close_polygons(grid.face_node_connectivity, grid.fill_value)
+    faces, _ = close_polygons(grid.face_node_connectivity)
     vertices = nodes[faces]
 
     # PolyCollection takes a norm, but not vmin, vmax.
diff --git a/xugrid/regrid/unstructured.py b/xugrid/regrid/unstructured.py
index ab8a6a3e6..b823af4ad 100644
--- a/xugrid/regrid/unstructured.py
+++ b/xugrid/regrid/unstructured.py
@@ -152,7 +152,6 @@ def barycentric(self, other):
             grid.centroids,
             edge_face_connectivity=grid.edge_face_connectivity,
             edge_node_connectivity=grid.edge_node_connectivity,
-            fill_value=grid.fill_value,
             add_exterior=True,
             add_vertices=True,
             skip_concave=True,
diff --git a/xugrid/ugrid/connectivity.py b/xugrid/ugrid/connectivity.py
index 82d66c38b..d5088554b 100644
--- a/xugrid/ugrid/connectivity.py
+++ b/xugrid/ugrid/connectivity.py
@@ -7,7 +7,14 @@
 import pandas as pd
 from scipy import sparse
 
-from xugrid.constants import BoolArray, FloatArray, IntArray, IntDType, SparseMatrix
+from xugrid.constants import (
+    FILL_VALUE,
+    BoolArray,
+    FloatArray,
+    IntArray,
+    IntDType,
+    SparseMatrix,
+)
 
 
 def argsort_rows(array: np.ndarray) -> IntArray:
@@ -215,10 +222,10 @@ def contract_vertices(A: sparse.csr_matrix, indices: IntArray) -> IntArray:
 
 # Conversion between dense and sparse
 # -----------------------------------
-def _to_ij(conn: IntArray, fill_value: int, invert: bool) -> Tuple[IntArray, IntArray]:
+def _to_ij(conn: IntArray, invert: bool) -> Tuple[IntArray, IntArray]:
     n, m = conn.shape
     j = conn.ravel()
-    valid = j != fill_value
+    valid = j != FILL_VALUE
     i = np.repeat(np.arange(n), m)[valid]
     j = j[valid]
     if invert:
@@ -227,10 +234,8 @@ def _to_ij(conn: IntArray, fill_value: int, invert: bool) -> Tuple[IntArray, Int
         return i, j
 
 
-def _to_sparse(
-    conn: IntArray, fill_value: int, invert: bool, sort_indices: bool
-) -> sparse.csr_matrix:
-    i, j = _to_ij(conn, fill_value, invert)
+def _to_sparse(conn: IntArray, invert: bool, sort_indices: bool) -> sparse.csr_matrix:
+    i, j = _to_ij(conn, invert)
     coo_content = (j, (i, j))
     coo_matrix = sparse.coo_matrix(coo_content)
     csr_matrix = coo_matrix.tocsr()
@@ -270,13 +275,11 @@ def ragged_index(n: int, m: int, m_per_row: IntArray) -> BoolArray:
     return (column_number.T < m_per_row).T
 
 
-def to_sparse(
-    conn: IntArray, fill_value: int, sort_indices: bool = True
-) -> sparse.csr_matrix:
-    return _to_sparse(conn, fill_value, invert=False, sort_indices=sort_indices)
+def to_sparse(conn: IntArray, sort_indices: bool = True) -> sparse.csr_matrix:
+    return _to_sparse(conn, invert=False, sort_indices=sort_indices)
 
 
-def to_dense(conn: SparseMatrix, fill_value: int, n_columns: int = None) -> IntArray:
+def to_dense(conn: SparseMatrix, n_columns: int = None) -> IntArray:
     n, _ = conn.shape
     m_per_row = conn.getnnz(axis=1)
     m = m_per_row.max()
@@ -296,7 +299,7 @@ def to_dense(conn: SparseMatrix, fill_value: int, n_columns: int = None) -> IntA
         valid = slice(None)  # a[:] equals a[slice(None)]
     else:
         valid = ragged_index(n, m, m_per_row).ravel()
-        flat_conn[~valid] = fill_value
+        flat_conn[~valid] = FILL_VALUE
 
     if isinstance(conn, sparse.csr_matrix):
         flat_conn[valid] = conn.indices
@@ -310,18 +313,14 @@ def to_dense(conn: SparseMatrix, fill_value: int, n_columns: int = None) -> IntA
 # Inverting connectivities
 # ------------------------
 def invert_dense_to_sparse(
-    conn: IntArray, fill_value: int, sort_indices: bool = True
+    conn: IntArray, sort_indices: bool = True
 ) -> sparse.csr_matrix:
-    return _to_sparse(conn, fill_value, invert=True, sort_indices=sort_indices)
+    return _to_sparse(conn, invert=True, sort_indices=sort_indices)
 
 
-def invert_dense(
-    conn: IntArray, fill_value: int, sort_indices: bool = True
-) -> IntArray:
-    sparse_inverted = _to_sparse(
-        conn, fill_value, invert=True, sort_indices=sort_indices
-    )
-    return to_dense(sparse_inverted, fill_value)
+def invert_dense(conn: IntArray, sort_indices: bool = True) -> IntArray:
+    sparse_inverted = _to_sparse(conn, invert=True, sort_indices=sort_indices)
+    return to_dense(sparse_inverted)
 
 
 def invert_sparse(conn: sparse.csr_matrix) -> sparse.csr_matrix:
@@ -333,9 +332,9 @@ def invert_sparse(conn: sparse.csr_matrix) -> sparse.csr_matrix:
     return inverted.tocsr()
 
 
-def invert_sparse_to_dense(conn: sparse.csr_matrix, fill_value: int) -> IntArray:
+def invert_sparse_to_dense(conn: sparse.csr_matrix) -> IntArray:
     inverted = invert_sparse(conn)
-    return to_dense(inverted, fill_value)
+    return to_dense(inverted)
 
 
 # Renumbering
@@ -353,51 +352,46 @@ def _renumber(a: IntArray) -> IntArray:
     return dense.reshape(a.shape)
 
 
-def renumber(a: IntArray, fill_value: int = None):
-    if fill_value is None:
-        return _renumber(a)
-
-    valid = a != fill_value
-    renumbered = np.full_like(a, fill_value)
+def renumber(a: IntArray):
+    valid = a != FILL_VALUE
+    renumbered = np.full_like(a, FILL_VALUE)
     renumbered[valid] = _renumber(a[valid])
     return renumbered
 
 
-def close_polygons(face_node_connectivity: IntArray, fill_value: int) -> IntArray:
+def close_polygons(face_node_connectivity: IntArray) -> IntArray:
     # Wrap around and create closed polygon: put the first node at the end of the row
     # In case of fill values, replace all fill values
     n, m = face_node_connectivity.shape
-    closed = np.full((n, m + 1), fill_value, dtype=IntDType)
+    closed = np.full((n, m + 1), FILL_VALUE, dtype=IntDType)
     closed[:, :-1] = face_node_connectivity
     first_node = face_node_connectivity[:, 0]
     # Identify fill value, and replace by first node also
-    isfill = closed == fill_value
+    isfill = closed == FILL_VALUE
     closed.ravel()[isfill.ravel()] = np.repeat(first_node, isfill.sum(axis=1))
     return closed, isfill
 
 
-def reverse_orientation(face_node_connectivity: IntArray, fill_value: int):
+def reverse_orientation(face_node_connectivity: IntArray):
     # We cannot simply reverse the rows with [:, ::-1], since there may be fill
     # values present.
     reversed_orientation = face_node_connectivity.copy()
     in_reverse = face_node_connectivity[:, ::-1]
-    in_reverse = in_reverse[in_reverse != fill_value]
-    replace = face_node_connectivity != fill_value
+    in_reverse = in_reverse[in_reverse != FILL_VALUE]
+    replace = face_node_connectivity != FILL_VALUE
     reversed_orientation[replace] = in_reverse
     return reversed_orientation
 
 
-def counterclockwise(
-    face_node_connectivity: IntArray, fill_value: int, nodes: FloatArray
-) -> IntArray:
+def counterclockwise(face_node_connectivity: IntArray, nodes: FloatArray) -> IntArray:
     # TODO: in case of "periodic grids", we need to wrap around the grid.
-    closed, _ = close_polygons(face_node_connectivity, fill_value)
+    closed, _ = close_polygons(face_node_connectivity)
     p = nodes[closed]
     dxy = np.diff(p, axis=1)
     reverse = (np.cross(dxy[:, :-1], dxy[:, 1:])).sum(axis=1) < 0
     ccw = face_node_connectivity.copy()
     if reverse.any():
-        ccw[reverse] = reverse_orientation(face_node_connectivity[reverse], fill_value)
+        ccw[reverse] = reverse_orientation(face_node_connectivity[reverse])
     return ccw
 
 
@@ -405,24 +399,22 @@ def counterclockwise(
 # ----------------------
 def boundary_node_connectivity(
     edge_face_connectivity: IntArray,
-    fill_value: int,
     edge_node_connectivity: IntArray,
 ) -> IntArray:
     """Is a subset of the edge_node_connectivity"""
-    is_boundary = (edge_face_connectivity == fill_value).any(axis=1)
+    is_boundary = (edge_face_connectivity == FILL_VALUE).any(axis=1)
     return edge_node_connectivity[is_boundary]
 
 
 def edge_connectivity(
     face_node_connectivity: IntArray,
-    fill_value: int,
     edge_node_connectivity=None,
 ) -> Tuple[IntArray, IntArray]:
     """Derive new edge_node_connectivity and face_edge_connectivity."""
     prior = edge_node_connectivity
     n, m = face_node_connectivity.shape
     # Close the polygons: [0 1 2 3] -> [0 1 2 3 0]
-    closed, isfill = close_polygons(face_node_connectivity, fill_value)
+    closed, isfill = close_polygons(face_node_connectivity)
     # Allocate array for edge_node_connectivity: includes duplicate edges
     edge_node_connectivity = np.empty((n * m, 2), dtype=IntDType)
     edge_node_connectivity[:, 0] = closed[:, :-1].ravel()
@@ -449,7 +441,7 @@ def edge_connectivity(
         edge_node_connectivity = prior
 
     # Create face_edge_connectivity
-    face_edge_connectivity = np.full((n, m), fill_value, dtype=np.int64)
+    face_edge_connectivity = np.full((n, m), FILL_VALUE, dtype=np.int64)
     isnode = ~isfill[:, :-1]
     face_edge_connectivity[isnode] = inverse_indices
     return edge_node_connectivity, face_edge_connectivity
@@ -457,14 +449,13 @@ def edge_connectivity(
 
 def validate_edge_node_connectivity(
     face_node_connectivity: IntArray,
-    fill_value: int,
     edge_node_connectivity: IntArray,
 ) -> BoolArray:
     """
     * Is the edge defined by the face_node_connectivity?
     * Are the edges unique?
     """
-    new, _ = edge_connectivity(face_node_connectivity, fill_value)
+    new, _ = edge_connectivity(face_node_connectivity)
     old = np.sort(edge_node_connectivity, axis=1)
 
     new1d = new.astype(np.int32).view(np.int64).ravel()
@@ -485,7 +476,6 @@ def validate_edge_node_connectivity(
 
 def face_face_connectivity(
     edge_face_connectivity: IntArray,
-    fill_value: int,
 ) -> sparse.csr_matrix:
     """
     Derive face to face connectivity.
@@ -495,7 +485,7 @@ def face_face_connectivity(
     """
     i = edge_face_connectivity[:, 0]
     j = edge_face_connectivity[:, 1]
-    is_connection = j != fill_value
+    is_connection = j != FILL_VALUE
     i = i[is_connection]
     j = j[is_connection]
     edge_index = np.arange(len(edge_face_connectivity))[is_connection]
@@ -557,12 +547,11 @@ def structured_connectivity(active: IntArray) -> AdjacencyMatrix:
 
 def perimeter(
     face_node_connectivity: IntArray,
-    fill_value: int,
     node_x: FloatArray,
     node_y: FloatArray,
 ):
     nodes = np.column_stack([node_x, node_y])
-    closed, _ = close_polygons(face_node_connectivity, fill_value)
+    closed, _ = close_polygons(face_node_connectivity)
     coordinates = nodes[closed]
     # Shift coordinates to avoid precision loss
     xy0 = coordinates[:, 0]
@@ -583,19 +572,17 @@ def area_from_coordinates(
 
 def area(
     face_node_connectivity: IntArray,
-    fill_value: int,
     node_x: FloatArray,
     node_y: FloatArray,
 ) -> FloatArray:
     nodes = np.column_stack([node_x, node_y])
-    closed, _ = close_polygons(face_node_connectivity, fill_value)
+    closed, _ = close_polygons(face_node_connectivity)
     coordinates = nodes[closed]
     return area_from_coordinates(coordinates)
 
 
 def centroids(
     face_node_connectivity: IntArray,
-    fill_value: int,
     node_x: FloatArray,
     node_y: FloatArray,
 ) -> FloatArray:
@@ -610,7 +597,7 @@ def centroids(
         # This is mathematically equivalent to triangulating, computing triangle centroids
         # and computing the area weighted average of those centroids
         centroid_coordinates = np.empty((n_face, 2), dtype=np.float64)
-        closed, _ = close_polygons(face_node_connectivity, fill_value)
+        closed, _ = close_polygons(face_node_connectivity)
         coordinates = nodes[closed]
         # Express coordinates relative to first node
         xy0 = coordinates[:, 0]
@@ -644,7 +631,6 @@ def _circumcenters_triangle(xxx: FloatArray, yyy: FloatArray):
 
 def circumcenters(
     face_node_connectivity: IntArray,
-    fill_value: int,
     node_x: FloatArray,
     node_y: FloatArray,
 ) -> FloatArray:
@@ -684,16 +670,14 @@ def _triangulate(i: IntArray, j: IntArray, n_triangle_per_row: IntArray) -> IntA
     return triangles
 
 
-def triangulate_dense(
-    face_node_connectivity: IntArray, fill_value: int
-) -> Tuple[IntArray, IntArray]:
+def triangulate_dense(face_node_connectivity: IntArray) -> Tuple[IntArray, IntArray]:
     n_face, n_max = face_node_connectivity.shape
 
     if n_max == 3:
         triangles = face_node_connectivity.copy()
         return triangles, np.arange(n_face)
 
-    valid = face_node_connectivity != fill_value
+    valid = face_node_connectivity != FILL_VALUE
     n_per_row = valid.sum(axis=1)
     n_triangle_per_row = n_per_row - 2
     i = np.repeat(np.arange(n_face), n_per_row)
@@ -727,9 +711,7 @@ def triangulate_coo(
     return triangles, triangle_face_connectivity
 
 
-def triangulate(
-    face_node_connectivity, fill_value: int | None = None
-) -> Tuple[IntArray, IntArray]:
+def triangulate(face_node_connectivity) -> Tuple[IntArray, IntArray]:
     """
     Convert polygons into its constituent triangles.
 
@@ -747,9 +729,7 @@ def triangulate(
     triangle_face_connectivity: ndarray of integers with shape ``(n_triangle,)``
     """
     if isinstance(face_node_connectivity, IntArray):
-        if fill_value is None:
-            raise ValueError("fill_value is required for dense connectivity")
-        return triangulate_dense(face_node_connectivity, fill_value)
+        return triangulate_dense(face_node_connectivity)
     elif isinstance(face_node_connectivity, sparse.coo_matrix):
         return triangulate_coo(face_node_connectivity)
     else:
diff --git a/xugrid/ugrid/partitioning.py b/xugrid/ugrid/partitioning.py
index e2a64e038..de6e253ae 100644
--- a/xugrid/ugrid/partitioning.py
+++ b/xugrid/ugrid/partitioning.py
@@ -6,7 +6,7 @@
 import numpy as np
 import xarray as xr
 
-from xugrid.constants import IntArray, IntDType
+from xugrid.constants import FILL_VALUE, IntArray, IntDType
 from xugrid.core.wrap import UgridDataArray, UgridDataset
 from xugrid.ugrid.connectivity import renumber
 from xugrid.ugrid.ugridbase import UgridType
@@ -115,17 +115,17 @@ def _merge_connectivity(gathered, slices):
     return merged, indexes
 
 
-def merge_faces(grids, node_inverse, fill_value: int = -1):
+def merge_faces(grids, node_inverse):
     node_offsets = tuple(accumulate([0] + [grid.n_node for grid in grids]))
     n_face = [grid.n_face for grid in grids]
     n_max_node = max(grid.n_max_node_per_face for grid in grids)
     slices = (0,) + tuple(accumulate(n_face))
 
-    all_faces = np.full((sum(n_face), n_max_node), fill_value, dtype=IntDType)
+    all_faces = np.full((sum(n_face), n_max_node), FILL_VALUE, dtype=IntDType)
     for grid, face_offset, node_offset in zip(grids, slices, node_offsets):
         faces = grid.face_node_connectivity
         n_face, n_node_per_face = faces.shape
-        valid = faces != grid.fill_value
+        valid = faces != FILL_VALUE
         all_faces[face_offset : face_offset + n_face, :n_node_per_face][
             valid
         ] = node_inverse[faces[valid] + node_offset]
diff --git a/xugrid/ugrid/polygonize.py b/xugrid/ugrid/polygonize.py
index 8a88b5a69..6e19c74a2 100644
--- a/xugrid/ugrid/polygonize.py
+++ b/xugrid/ugrid/polygonize.py
@@ -3,7 +3,7 @@
 import numpy as np
 from scipy import sparse
 
-from xugrid.constants import IntArray
+from xugrid.constants import FILL_VALUE, IntArray
 
 
 def _bbox_area(bounds):
@@ -11,7 +11,7 @@ def _bbox_area(bounds):
 
 
 def _classify(
-    fill_value: int, i: IntArray, j: IntArray, face_values: np.ndarray
+    i: IntArray, j: IntArray, face_values: np.ndarray
 ) -> Tuple[int, IntArray]:
     """
     Find out how many discrete polygons are created. Identify the connectivity,
@@ -19,8 +19,6 @@ def _classify(
 
     Parameters
     ----------
-    fill_value: int
-        Fill value in j: marks exterior edges.
     i: np.ndarray of int
         First face of the edge.
     j: np.ndarray of int
@@ -40,7 +38,7 @@ def _classify(
     # For labelling, only those parts of the mesh that have the same value
     # should be connected with each other.
     # Since we dropped NaN values before, we needn't worry about those.
-    is_connection = (i != fill_value) & (j != fill_value) & (vi == vj)
+    is_connection = (i != FILL_VALUE) & (j != FILL_VALUE) & (vi == vj)
     i = i[is_connection]
     j = j[is_connection]
     ij = np.concatenate([i, j])
@@ -94,8 +92,7 @@ def polygonize(uda: "UgridDataArray") -> "gpd.GeoDataFrame":  # type: ignore # n
     face_values = dropped.to_numpy()
     grid = dropped.ugrid.grid
     i, j = grid.edge_face_connectivity.T
-    fill_value = grid.fill_value
-    n_polygon, polygon_id = _classify(fill_value, i, j, face_values)
+    n_polygon, polygon_id = _classify(i, j, face_values)
 
     # Now we identify for each label the subset of edges. These are the
     # "exterior" edges: either the exterior edge of the mesh identified by a
@@ -106,8 +103,8 @@ def polygonize(uda: "UgridDataArray") -> "gpd.GeoDataFrame":  # type: ignore # n
     vj = polygon_id[j]
     # Ensure that no result thas has been created by indexing with the
     # fill_value remains. Since polygon_id starts counting a 0, we may use -1.
-    vi[i == fill_value] = -1
-    vj[j == fill_value] = -1
+    vi[i == FILL_VALUE] = FILL_VALUE
+    vj[j == FILL_VALUE] = FILL_VALUE
     boundary = vi != vj
 
     polygons = []
diff --git a/xugrid/ugrid/snapping.py b/xugrid/ugrid/snapping.py
index 53b93e29e..b8a842144 100644
--- a/xugrid/ugrid/snapping.py
+++ b/xugrid/ugrid/snapping.py
@@ -402,7 +402,7 @@ def create_snap_to_grid_dataframe(
     vertices = topology.node_coordinates
     edge_centroids = topology.edge_coordinates
     face_edge_connectivity = topology.face_edge_connectivity
-    A = connectivity.to_sparse(face_edge_connectivity, fill_value=-1)
+    A = connectivity.to_sparse(face_edge_connectivity)
     n, m = A.shape
     face_edge_connectivity = AdjacencyMatrix(A.indices, A.indptr, A.nnz, n, m)
 
diff --git a/xugrid/ugrid/ugrid1d.py b/xugrid/ugrid/ugrid1d.py
index 16ab0dcd6..4d8f74bf9 100644
--- a/xugrid/ugrid/ugrid1d.py
+++ b/xugrid/ugrid/ugrid1d.py
@@ -9,6 +9,7 @@
 import xugrid
 from xugrid import conversion
 from xugrid.constants import (
+    FILL_VALUE,
     BoolArray,
     FloatArray,
     FloatDType,
@@ -133,7 +134,6 @@ def from_dataset(cls, dataset: xr.Dataset, topology: str = None):
             + list(connectivity.values())
             + list(chain.from_iterable(chain.from_iterable(coordinates.values())))
         )
-        fill_value = -1
 
         # Take the first coordinates by default.
         # They can be reset with .set_node_coords()
@@ -144,7 +144,7 @@ def from_dataset(cls, dataset: xr.Dataset, topology: str = None):
 
         edge_nodes = connectivity["edge_node_connectivity"]
         edge_node_connectivity = cls._prepare_connectivity(
-            ds[edge_nodes], fill_value, dtype=IntDType
+            ds[edge_nodes], FILL_VALUE, dtype=IntDType
         ).to_numpy()
 
         indexes["node_x"] = x_index
@@ -154,7 +154,7 @@ def from_dataset(cls, dataset: xr.Dataset, topology: str = None):
         return cls(
             node_x_coordinates,
             node_y_coordinates,
-            fill_value,
+            FILL_VALUE,
             edge_node_connectivity,
             name=topology,
             dataset=dataset[ugrid_vars],
@@ -211,7 +211,7 @@ def from_meshkernel(
         return cls(
             mesh.node_x,
             mesh.node_y,
-            fill_value=-1,
+            fill_value=FILL_VALUE,
             edge_node_connectivity=mesh.edge_nodes.reshape((-1, 2)),
             name=name,
             projected=projected,
@@ -389,8 +389,7 @@ def from_shapely(geometry: LineArray, crs=None) -> "Ugrid1d":
             )
 
         x, y, edge_node_connectivity = conversion.linestrings_to_edges(geometry)
-        fill_value = -1
-        return Ugrid1d(x, y, fill_value, edge_node_connectivity, crs=crs)
+        return Ugrid1d(x, y, FILL_VALUE, edge_node_connectivity, crs=crs)
 
     def to_pygeos(self, dim):
         from warnings import warn
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index 1ea8c936c..377c1b3e9 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -16,6 +16,7 @@
 from xugrid import conversion
 from xugrid import meshkernel_utils as mku
 from xugrid.constants import (
+    FILL_VALUE,
     BoolArray,
     FloatArray,
     FloatDType,
@@ -104,17 +105,18 @@ def __init__(
         self.projected = projected
 
         if isinstance(face_node_connectivity, np.ndarray):
-            face_node_connectivity = face_node_connectivity
+            self.face_node_connectivity = face_node_connectivity.copy()
         elif isinstance(face_node_connectivity, (coo_matrix, csr_matrix)):
-            face_node_connectivity = connectivity.to_dense(
-                face_node_connectivity, fill_value
-            )
+            self.face_node_connectivity = connectivity.to_dense(face_node_connectivity)
         else:
             raise TypeError(
                 "face_node_connectivity should be an array of integers or a sparse matrix"
             )
 
-        self.face_node_connectivity = face_node_connectivity
+        # Ensure the fill value is FILL_VALUE (-1)
+        self.face_node_connectivity[
+            self.face_node_connectivity == self.fill_value
+        ] = FILL_VALUE
 
         # TODO: do this in validation instead. While UGRID conventions demand it,
         # where does it go wrong?
@@ -225,15 +227,14 @@ def from_meshkernel(
         """
         n_face = len(mesh.nodes_per_face)
         n_max_node = mesh.nodes_per_face.max()
-        fill_value = -1
-        face_node_connectivity = np.full((n_face, n_max_node), fill_value)
+        face_node_connectivity = np.full((n_face, n_max_node), FILL_VALUE)
         isnode = connectivity.ragged_index(n_face, n_max_node, mesh.nodes_per_face)
         face_node_connectivity[isnode] = mesh.face_nodes
         edge_node_connectivity = np.reshape(mesh.edge_nodes, (-1, 2))
         return cls(
             node_x=mesh.node_x,
             node_y=mesh.node_y,
-            fill_value=fill_value,
+            fill_value=FILL_VALUE,
             face_node_connectivity=face_node_connectivity,
             edge_node_connectivity=edge_node_connectivity,
             name=name,
@@ -330,14 +331,14 @@ def to_dataset(
         data_vars = {
             self.name: 0,
             face_nodes: xr.DataArray(
-                data=self.face_node_connectivity,
+                data=self._set_fillvalue(self.face_node_connectivity),
                 attrs=face_nodes_attrs,
                 dims=(self.face_dimension, nmax_node_dim),
             ),
         }
         if self.edge_node_connectivity is not None or optional_attributes:
             data_vars[edge_nodes] = xr.DataArray(
-                data=self.edge_node_connectivity,
+                data=self.edge_node_connectivity,  # has no fill values
                 attrs=edge_nodes_attrs,
                 dims=(self.edge_dimension, "two"),
             )
@@ -346,28 +347,29 @@ def to_dataset(
             face_faces, face_faces_attrs = self._get_name_and_attrs("face_face")
             edge_faces, edge_faces_attrs = self._get_name_and_attrs("edge_face")
             bound_nodes, bound_nodes_attrs = self._get_name_and_attrs("boundary_node")
-            fill_value = self.fill_value
             boundary_edge_dim = self._attrs["boundary_edge_dimension"]
 
             data_vars[face_edges] = xr.DataArray(
-                data=self.face_edge_connectivity,
+                data=self._set_fillvalue(self.face_edge_connectivity),
                 attrs=face_edges_attrs,
                 dims=(self.face_dimension, nmax_node_dim),
             )
             data_vars[face_faces] = xr.DataArray(
-                data=connectivity.to_dense(
-                    self.face_face_connectivity, fill_value, self.n_max_node_per_face
+                data=self._set_fillvalue(
+                    connectivity.to_dense(
+                        self.face_face_connectivity, self.n_max_node_per_face
+                    )
                 ),
                 attrs=face_faces_attrs,
                 dims=(self.face_dimension, nmax_node_dim),
             )
             data_vars[edge_faces] = xr.DataArray(
-                data=self.edge_face_connectivity,
+                data=self._set_fillvalue(self.edge_face_connectivity),
                 attrs=edge_faces_attrs,
                 dims=(self.edge_dimension, "two"),
             )
             data_vars[bound_nodes] = xr.DataArray(
-                data=self.boundary_node_connectivity,
+                data=self._set_fillvalue(self.boundary_node_connectivity),
                 attrs=bound_nodes_attrs,
                 dims=(boundary_edge_dim, "two"),
             )
@@ -407,7 +409,7 @@ def n_max_node_per_face(self) -> int:
 
     @property
     def n_node_per_face(self) -> IntArray:
-        return (self.face_node_connectivity != self.fill_value).sum(axis=1)
+        return (self.face_node_connectivity != FILL_VALUE).sum(axis=1)
 
     @property
     def core_dimension(self):
@@ -460,7 +462,6 @@ def _edge_connectivity(self):
             self._face_edge_connectivity,
         ) = connectivity.edge_connectivity(
             self.face_node_connectivity,
-            self.fill_value,
             self._edge_node_connectivity,
         )
 
@@ -507,7 +508,6 @@ def boundary_node_connectivity(self) -> IntArray:
         if self._boundary_node_connectivity is None:
             self._boundary_node_connectivity = connectivity.boundary_node_connectivity(
                 self.edge_face_connectivity,
-                self.fill_value,
                 self.edge_node_connectivity,
             )
         return self._boundary_node_connectivity
@@ -524,7 +524,6 @@ def centroids(self) -> FloatArray:
         if self._centroids is None:
             self._centroids = connectivity.centroids(
                 self.face_node_connectivity,
-                self.fill_value,
                 self.node_x,
                 self.node_y,
             )
@@ -539,7 +538,6 @@ def circumcenters(self):
         if self._circumcenters is None:
             self._circumcenters = connectivity.circumcenters(
                 self.face_node_connectivity,
-                self.fill_value,
                 self.node_x,
                 self.node_y,
             )
@@ -551,7 +549,6 @@ def area(self) -> FloatArray:
         if self._area is None:
             self._area = connectivity.area(
                 self.face_node_connectivity,
-                self.fill_value,
                 self.node_x,
                 self.node_y,
             )
@@ -563,7 +560,6 @@ def perimeter(self) -> FloatArray:
         if self._perimeter is None:
             self._perimeter = connectivity.perimeter(
                 self.face_node_connectivity,
-                self.fill_value,
                 self.node_x,
                 self.node_y,
             )
@@ -581,7 +577,7 @@ def face_bounds(self):
         """
         x = self.node_x[self.face_node_connectivity]
         y = self.node_y[self.face_node_connectivity]
-        isfill = self.face_node_connectivity == self.fill_value
+        isfill = self.face_node_connectivity == FILL_VALUE
         x[isfill] = np.nan
         y[isfill] = np.nan
         return np.column_stack(
@@ -628,7 +624,7 @@ def face_node_coordinates(self) -> FloatArray:
         coords = np.full(
             (self.n_face, self.n_max_node_per_face, 2), np.nan, dtype=FloatDType
         )
-        is_node = self.face_node_connectivity != self.fill_value
+        is_node = self.face_node_connectivity != FILL_VALUE
         index = self.face_node_connectivity[is_node]
         coords[is_node, :] = self.node_coordinates[index]
         return coords
@@ -639,7 +635,7 @@ def edge_face_connectivity(self) -> IntArray:
         Edge to face connectivity. An edge may belong to a single face
         (exterior edge), or it may be shared by two faces (interior edge).
 
-        An exterior edge will contain a ``fill_value`` for the second column.
+        An exterior edge will contain a FILL_VALUE of -1 for the second column.
 
         Returns
         -------
@@ -647,7 +643,7 @@ def edge_face_connectivity(self) -> IntArray:
         """
         if self._edge_face_connectivity is None:
             self._edge_face_connectivity = connectivity.invert_dense(
-                self.face_edge_connectivity, self.fill_value
+                self.face_edge_connectivity
             )
         return self._edge_face_connectivity
 
@@ -667,7 +663,7 @@ def face_face_connectivity(self) -> csr_matrix:
         """
         if self._face_face_connectivity is None:
             self._face_face_connectivity = connectivity.face_face_connectivity(
-                self.edge_face_connectivity, self.fill_value
+                self.edge_face_connectivity
             )
         return self._face_face_connectivity
 
@@ -682,7 +678,7 @@ def node_face_connectivity(self):
         """
         if self._node_face_connectivity is None:
             self._node_face_connectivity = connectivity.invert_dense_to_sparse(
-                self.face_node_connectivity, self.fill_value
+                self.face_node_connectivity
             )
         return self._node_face_connectivity
 
@@ -739,7 +735,7 @@ def mesh(self) -> "mk.Mesh2d":  # type: ignore # noqa
         import meshkernel as mk
 
         edge_nodes = self.edge_node_connectivity.ravel().astype(np.int32)
-        is_node = self.face_node_connectivity != self.fill_value
+        is_node = self.face_node_connectivity != FILL_VALUE
         nodes_per_face = is_node.sum(axis=1).astype(np.int32)
         face_nodes = self.face_node_connectivity[is_node].ravel().astype(np.int32)
 
@@ -792,7 +788,6 @@ def voronoi_topology(self):
                 self.centroids,
                 self.edge_face_connectivity,
                 self.edge_node_connectivity,
-                self.fill_value,
                 add_exterior=True,
                 add_vertices=False,
             )
@@ -816,7 +811,7 @@ def centroid_triangulation(self):
         """
         if self._centroid_triangulation is None:
             nodes, faces, face_index = self.voronoi_topology
-            triangles, _ = connectivity.triangulate(faces, self.fill_value)
+            triangles, _ = connectivity.triangulate(faces)
             triangulation = (nodes[:, 0].copy(), nodes[:, 1].copy(), triangles)
             self._centroid_triangulation = (triangulation, face_index)
         return self._centroid_triangulation
@@ -835,7 +830,7 @@ def triangulation(self):
         """
         if self._triangulation is None:
             triangles, triangle_face_connectivity = connectivity.triangulate(
-                self.face_node_connectivity, self.fill_value
+                self.face_node_connectivity
             )
             triangulation = (self.node_x, self.node_y, triangles)
             self._triangulation = (triangulation, triangle_face_connectivity)
@@ -851,7 +846,7 @@ def exterior_edges(self) -> IntArray:
         edge_index: 1d array of integers
         """
         # Numpy argwhere doesn't return a 1D array
-        return np.nonzero(self.edge_face_connectivity[:, 1] == self.fill_value)[0]
+        return np.nonzero(self.edge_face_connectivity[:, 1] == FILL_VALUE)[0]
 
     @property
     def exterior_faces(self) -> IntArray:
@@ -864,7 +859,7 @@ def exterior_faces(self) -> IntArray:
         """
         exterior_edges = self.exterior_edges
         exterior_faces = self.edge_face_connectivity[exterior_edges].ravel()
-        return np.unique(exterior_faces[exterior_faces != self.fill_value])
+        return np.unique(exterior_faces[exterior_faces != FILL_VALUE])
 
     @property
     def celltree(self):
@@ -875,7 +870,7 @@ def celltree(self):
         """
         if self._celltree is None:
             self._celltree = CellTree2d(
-                self.node_coordinates, self.face_node_connectivity, self.fill_value
+                self.node_coordinates, self.face_node_connectivity, FILL_VALUE
             )
         return self._celltree
 
@@ -904,7 +899,6 @@ def validate_edge_node_connectivity(self):
         """
         return connectivity.validate_edge_node_connectivity(
             self.face_node_connectivity,
-            self.fill_value,
             self.edge_node_connectivity,
         )
 
@@ -1121,8 +1115,8 @@ def topology_subset(
         index = face_index.to_numpy()
         face_subset = self.face_node_connectivity[index]
         node_index = np.unique(face_subset.ravel())
-        node_index = node_index[node_index != self.fill_value]
-        new_faces = connectivity.renumber(face_subset, self.fill_value)
+        node_index = node_index[node_index != FILL_VALUE]
+        new_faces = connectivity.renumber(face_subset)
         node_x = self.node_x[node_index]
         node_y = self.node_y[node_index]
 
@@ -1130,7 +1124,7 @@ def topology_subset(
         new_edges = None
         if self.edge_node_connectivity is not None:
             edge_index = np.unique(self.face_edge_connectivity[index].ravel())
-            edge_index = edge_index[edge_index != self.fill_value]
+            edge_index = edge_index[edge_index != FILL_VALUE]
             edge_subset = self.edge_node_connectivity[edge_index]
             new_edges = connectivity.renumber(edge_subset)
 
@@ -1214,7 +1208,7 @@ def isel(self, indexers=None, return_index=False, **indexers_kwargs):
         if edgedim in indexers:
             edge_index = indexers[edgedim]
             index = np.unique(self.edge_face_connectivity[edge_index])
-            face_index[edgedim] = index[index != self.fill_value]
+            face_index[edgedim] = index[index != FILL_VALUE]
         if facedim in indexers:
             face_index[facedim] = indexers[facedim]
 
@@ -1604,11 +1598,8 @@ def merge_partitions(
 
         # Grab a sample grid
         grid = next(iter(grids))
-        fill_value = grid.fill_value
         node_coordinates, node_indexes, node_inverse = partitioning.merge_nodes(grids)
-        new_faces, face_indexes = partitioning.merge_faces(
-            grids, node_inverse, fill_value
-        )
+        new_faces, face_indexes = partitioning.merge_faces(grids, node_inverse)
         indexes = {
             grid.node_dimension: node_indexes,
             grid.face_dimension: face_indexes,
@@ -1622,7 +1613,7 @@ def merge_partitions(
 
         merged_grid = Ugrid2d(
             *node_coordinates.T,
-            fill_value,
+            FILL_VALUE,
             new_faces,
             name=grid.name,
             edge_node_connectivity=new_edges,
@@ -1862,10 +1853,8 @@ def triangulate(self):
         -------
         triangles: Ugrid2d
         """
-        triangles, _ = connectivity.triangulate(
-            self.face_node_connectivity, self.fill_value
-        )
-        return Ugrid2d(self.node_x, self.node_y, self.fill_value, triangles)
+        triangles, _ = connectivity.triangulate(self.face_node_connectivity)
+        return Ugrid2d(self.node_x, self.node_y, FILL_VALUE, triangles)
 
     def _tesselate_voronoi(self, centroids, add_exterior, add_vertices, skip_concave):
         if add_exterior:
@@ -1881,12 +1870,11 @@ def _tesselate_voronoi(self, centroids, add_exterior, add_vertices, skip_concave
             centroids,
             edge_face_connectivity,
             edge_node_connectivity,
-            self.fill_value,
             add_exterior,
             add_vertices,
             skip_concave,
         )
-        return Ugrid2d(vertices[:, 0], vertices[:, 1], self.fill_value, faces)
+        return Ugrid2d(vertices[:, 0], vertices[:, 1], FILL_VALUE, faces)
 
     def tesselate_centroidal_voronoi(
         self, add_exterior=True, add_vertices=True, skip_concave=False
@@ -2092,10 +2080,8 @@ def from_shapely(geometry: PolygonArray, crs=None) -> "Ugrid2d":
                 "geometry contains other types of geometries."
             )
 
-        x, y, face_node_connectivity, fill_value = conversion.polygons_to_faces(
-            geometry
-        )
-        return Ugrid2d(x, y, fill_value, face_node_connectivity, crs=crs)
+        x, y, face_node_connectivity = conversion.polygons_to_faces(geometry)
+        return Ugrid2d(x, y, FILL_VALUE, face_node_connectivity, crs=crs)
 
     @staticmethod
     def _from_intervals_helper(
@@ -2295,7 +2281,6 @@ def to_shapely(self, dim):
                 self.node_x,
                 self.node_y,
                 self.face_node_connectivity,
-                self.fill_value,
             )
         elif dim == self.node_dimension:
             return conversion.nodes_to_points(
diff --git a/xugrid/ugrid/ugridbase.py b/xugrid/ugrid/ugridbase.py
index 8a7295e6a..77c940f8f 100644
--- a/xugrid/ugrid/ugridbase.py
+++ b/xugrid/ugrid/ugridbase.py
@@ -10,7 +10,7 @@
 from numpy.typing import ArrayLike
 from scipy.sparse import csr_matrix
 
-from xugrid.constants import BoolArray, FloatArray, IntArray
+from xugrid.constants import FILL_VALUE, BoolArray, FloatArray, IntArray
 from xugrid.ugrid import connectivity, conventions
 
 
@@ -483,6 +483,12 @@ def _prepare_connectivity(
             raise ValueError("connectivity contains negative values")
         return da.copy(data=cast)
 
+    def _set_fillvalue(self, connectivity: IntArray) -> IntArray:
+        c = connectivity.copy()
+        if self.fill_value != FILL_VALUE:
+            c[c == FILL_VALUE] = self.fill_value
+        return c
+
     def _precheck(self, multi_index):
         dim, index = multi_index.popitem()
         for check_dim, check_index in multi_index.items():
@@ -646,7 +652,7 @@ def node_edge_connectivity(self) -> csr_matrix:
         """
         if self._node_edge_connectivity is None:
             self._node_edge_connectivity = connectivity.invert_dense_to_sparse(
-                self.edge_node_connectivity, self.fill_value
+                self.edge_node_connectivity
             )
         return self._node_edge_connectivity
 
diff --git a/xugrid/ugrid/voronoi.py b/xugrid/ugrid/voronoi.py
index 84dee3d7d..7edb057b9 100644
--- a/xugrid/ugrid/voronoi.py
+++ b/xugrid/ugrid/voronoi.py
@@ -20,7 +20,7 @@
 import pandas as pd
 from scipy import sparse
 
-from xugrid.constants import X_EPSILON, FloatArray, IntArray
+from xugrid.constants import FILL_VALUE, X_EPSILON, FloatArray, IntArray
 from xugrid.ugrid.connectivity import (
     area_from_coordinates,
     close_polygons,
@@ -61,7 +61,7 @@ def _create_face_node_connectivity(i: IntArray, j: IntArray) -> IntArray:
     n = len(n_vertex)
     m = n_vertex.max()
     index = ragged_index(n, m, n_vertex)
-    face_node_connectivity = np.full((n, m), -1)
+    face_node_connectivity = np.full((n, m), FILL_VALUE)
     face_node_connectivity[index] = j
     return face_node_connectivity
 
@@ -79,10 +79,9 @@ def interior_centroids(
     node_face_connectivity: sparse.csr_matrix,
     edge_face_connectivity: IntArray,
     edge_node_connectivity: IntArray,
-    fill_value: int,
 ):
     # Find exterior nodes associated with interior edges
-    is_exterior = edge_face_connectivity[:, 1] == fill_value
+    is_exterior = edge_face_connectivity[:, 1] == FILL_VALUE
     exterior_nodes = np.unique(edge_node_connectivity[is_exterior].ravel())
     m_per_node = node_face_connectivity.getnnz(axis=1)
     is_interior_only = m_per_node > 1
@@ -146,12 +145,11 @@ def _interpolate_between_projections(
 def exterior_vertices(
     edge_face_connectivity: IntArray,
     edge_node_connectivity: IntArray,
-    fill_value: int,
     vertices: FloatArray,
     centroids: FloatArray,
     add_vertices: bool,
 ):
-    is_exterior = edge_face_connectivity[:, 1] == fill_value
+    is_exterior = edge_face_connectivity[:, 1] == FILL_VALUE
     exterior_nodes = edge_node_connectivity[is_exterior]
     # For every exterior node, project the centroids to exterior edges
     edge_vertices = vertices[exterior_nodes]
@@ -211,7 +209,7 @@ def choose_convex(
     faces = np.full((n, m), -1)
     faces[index] = j
     # Close the polygons so we can easily compute areas.
-    closed, _ = close_polygons(faces, -1)
+    closed, _ = close_polygons(faces)
     # Make a copy and insert the original vertices.
     modified_nodes = nodes.copy()
     modified_nodes[-n_interpolated:] = original_vertices
@@ -232,7 +230,6 @@ def exterior_topology(
     edge_face_connectivity: IntArray,
     edge_node_connectivity: IntArray,
     node_face_connectivity: sparse.csr_matrix,
-    fill_value: int,
     vertices: FloatArray,
     centroids: FloatArray,
     add_vertices: bool,
@@ -278,7 +275,6 @@ def exterior_topology(
         node_face_connectivity,
         edge_face_connectivity,
         edge_node_connectivity,
-        fill_value,
     )
     i1, j1 = exterior_centroids(node_face_connectivity)
     (
@@ -291,7 +287,6 @@ def exterior_topology(
     ) = exterior_vertices(
         edge_face_connectivity,
         edge_node_connectivity,
-        fill_value,
         vertices,
         centroids,
         add_vertices,
@@ -337,7 +332,6 @@ def voronoi_topology(
     centroids: FloatArray,
     edge_face_connectivity: IntArray = None,
     edge_node_connectivity: IntArray = None,
-    fill_value: int = None,
     add_exterior: bool = False,
     add_vertices: bool = False,
     skip_concave: bool = False,
@@ -377,8 +371,6 @@ def voronoi_topology(
     centroids: ndarray of floats with shape ``(n_centroid, 2)``
     edge_face_connectivity: ndarray of integers with shape ``(n_edge, 2)``, optional
     edge_node_connectivity: ndarray of integers with shape ``(n_edge, 2)``, optional
-    fill_value: int, optional
-        Fill value for edge_face_connectivity.
     add_exterior: bool, optional
         Whether to consider exterior edges of the original mesh, or to consider
         exclusively centroids.
@@ -401,12 +393,9 @@ def voronoi_topology(
         interpolated.
     """
     if add_exterior:
-        if any(
-            arg is None
-            for arg in [edge_face_connectivity, edge_node_connectivity, fill_value]
-        ):
+        if any(arg is None for arg in [edge_face_connectivity, edge_node_connectivity]):
             raise ValueError(
-                "edge_face_connectivity, edge_node_connectivity, fill_value "
+                "edge_face_connectivity, edge_node_connectivity "
                 "must be provided if add_exterior is True."
             )
 
@@ -415,7 +404,7 @@ def voronoi_topology(
     # take any valid internal polygon we can construct: at least a triangle.
     ncol_per_row = node_face_connectivity.getnnz(axis=1)
     if add_exterior:
-        is_exterior = edge_face_connectivity[:, 1] == fill_value
+        is_exterior = edge_face_connectivity[:, 1] == FILL_VALUE
         exterior_nodes = edge_node_connectivity[is_exterior]
         valid = np.full(len(vertices), True)
         valid[exterior_nodes.ravel()] = False
@@ -449,7 +438,6 @@ def voronoi_topology(
             edge_face_connectivity,
             edge_node_connectivity,
             node_face_connectivity,
-            fill_value,
             vertices,
             centroids,
             add_vertices,