diff --git a/CHANGELOG.md b/CHANGELOG.md
index 289feb9890..3ec7700e5d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@
 
 ## Bug fixes
 
+- 3D unstructured-mesh point-in-domain tests no longer misclassify points lying exactly on a domain face, edge, or corner as outside, so interpolating or plotting a field at the domain boundary no longer blanks it to NaN. ([#5704](https://github.com/pybamm-team/PyBaMM/pull/5704))
 - `ElectrodeSOHSolver` now passes model options through, so hysteresis OCP branches are used. ([#5701](https://github.com/pybamm-team/PyBaMM/pull/5701))
 - Fixed a memory leak in `ElectrodeSOHSolver.theoretical_energy_integral`, which cached a new expression tree per call. ([#5695](https://github.com/pybamm-team/PyBaMM/pull/5695))
 - `BatchStudy.solve` no longer ignores its `solver` argument: previously the loop over study inputs shadowed it, so a caller-supplied solver was silently dropped. A solver from `BatchStudy(solvers=...)` still takes precedence. ([#5677](https://github.com/pybamm-team/PyBaMM/pull/5677))
diff --git a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
index f2d344fd44..37e3b66d64 100644
--- a/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
+++ b/packages/pybamm/src/pybamm/meshes/unstructured_submesh.py
@@ -681,7 +681,11 @@ def contains_points_3d(self, query_pts):
 
             winding[i] = 2.0 * np.arctan2(num, den).sum()
 
-        return winding > 2.0 * np.pi
+        # Interior points have winding 4*pi and exterior 0; a point lying on a
+        # face/edge/corner has a partial value (2*pi / pi / pi-over-2). Use a
+        # small positive threshold so boundary points count as inside rather
+        # than being masked to NaN, while true exterior (winding ~ 0) stays out.
+        return winding > 0.1
 
 
 # ======================================================================
diff --git a/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py b/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py
index f37772781d..417f1a4bd4 100644
--- a/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py
+++ b/packages/pybamm/tests/unit/test_meshes/test_unstructured_submesh.py
@@ -406,6 +406,25 @@ def test_contains_points_3d_hex_mesh(self):
         assert mesh.contains_points_3d(np.array([[0.5, 0.5, 0.5]]))[0]
         assert not mesh.contains_points_3d(np.array([[2.0, 2.0, 2.0]]))[0]
 
+    def test_contains_points_3d_boundary_points_are_inside(self):
+        """Points lying exactly on a face, edge, or corner of the domain must
+        count as inside, not be masked out (their winding number is a partial
+        2*pi / pi / pi-over-2 rather than the full 4*pi)."""
+        mesh = UnstructuredSubMesh(*_hex_grid(*[np.linspace(0, 1, 3)] * 3))
+        mesh.detect_box_boundaries()
+
+        on_or_inside = np.array(
+            [
+                [0.0, 0.0, 0.0],  # corner
+                [0.5, 0.0, 0.0],  # edge midpoint
+                [0.5, 0.5, 0.0],  # face centre
+                [0.5, 0.5, 0.5],  # interior
+            ]
+        )
+        outside = np.array([[2.0, 2.0, 2.0]])
+        assert mesh.contains_points_3d(on_or_inside).all()
+        assert not mesh.contains_points_3d(outside)[0]
+
     def test_optimize_ordering_single_cell_noop(self):
         """optimize_ordering with 1 cell returns without permuting."""
         nodes = np.array([[0, 0], [1, 0], [0, 1]], dtype=float)