diff -r a4137aedcb3a -r c3a4f4bb87f7 src/dolfinx_access/function.rs
--- a/src/dolfinx_access/function.rs	Thu Feb 26 09:32:12 2026 -0500
+++ b/src/dolfinx_access/function.rs	Wed Apr 22 22:32:00 2026 -0500
@@ -1,5 +1,5 @@
 /*!
-Python and C++ wrapper for Dolfinx Function<f64>
+Python and C++ wrapper for Dolfinx `Function<f64>`.
 */
 
 use super::ffi;
@@ -43,15 +43,21 @@
     cxx: *mut Function_f64,
 }
 
+/// References to Python side helper routines
 struct Helpers {
+    /// Dot product
     dot: Py<PyFunction>,
+    /// Squared 2-norm
     norm2_squared: Py<PyFunction>,
+    /// Squared distance
     dist2_squared: Py<PyFunction>,
 }
 
 use std::sync::OnceLock;
+/// A populated-once structure for referencing Python-side helper routines.
 static HELPERS: OnceLock<Helpers> = OnceLock::new();
 
+/// Gets [`Helpers`], populates [`HELPERS`] if needed.
 fn get_helpers(py: Python<'_>) -> PyResult<&Helpers> {
     HELPERS.get_or_try_init(|| {
         let extras = PyModule::import(py, "pointsource_pde.dolfinx_extras")?;
@@ -72,6 +78,7 @@
                 unsafe { ffi::drop_Function_f64(self.cxx) };
             }
         }
+        self.cxx = std::ptr::null_mut();
     }
 }
 
@@ -80,7 +87,7 @@
         DolfinxPyFunction_f64 {
             u: None, // It's a different new function, so no Python instance
             function_space: self.function_space.clone(),
-            owned: self.owned,
+            owned: true,
             cxx: unsafe { ffi::clone_Function_f64(self.cxx) },
         }
     }
@@ -135,7 +142,7 @@
 
     /// Returns a (possibly new) Python presentation, without consuming current object.
     ///
-    /// Please use this functionality through [`pyo3::types::IntoPyObject`].
+    /// Please use this functionality through [`pyo3::IntoPyObject`].
     pub(crate) fn _into_py_ref<'a>(
         &'a mut self,
         py: Python<'py>,
@@ -150,7 +157,7 @@
 
     /// Returns a (possibly new) Python presentation, consuming current object.
     ///
-    /// Please use this functionality through [`pyo3::types::IntoPyObject`].
+    /// Please use this functionality through [`pyo3::IntoPyObject`].
     pub(crate) fn _into_py(mut self, py: Python<'py>) -> PyResult<pyo3::Bound<'py, PyAny>> {
         self._make_py(py);
         self.cxx = std::ptr::null_mut();
@@ -181,6 +188,9 @@
     }
 }
 
+/// An extremely durty hack to create a Python-side Dolfinx `Function` referencing
+/// a C++-side Dolfinx `Function` that we already have. For some unfathomable reason,
+/// Dolfinx provies no way to do this cleanly.
 pub(crate) static CREATE_FUNCTION_HACK: &CStr = c_str!(
     "\
 from dolfinx.fem import Function
@@ -198,9 +208,12 @@
 "
 );
 
+/// Trait for points that can be expanded into Dolfinx' always-3D presentation
 trait ExpandCoordDolphinx {
+    /// Convert to always-3D presentation
     fn to_dolphinx(self) -> [f64; 3];
     #[allow(dead_code)]
+    /// Convert from always-3D presentation
     fn from_dolphinx(val: [f64; 3]) -> Self;
 }
 
@@ -243,6 +256,7 @@
     }
 }
 
+/// Implements [`DifferentiableImpl`]
 macro_rules! impl_diff {
     ($n:literal ; $($o:literal)+) => { $(
         impl<'py> DifferentiableImpl<Loc<$n, f64>> for DolfinxPyFunction_f64<'py, $n, $o, 1> {
@@ -274,6 +288,7 @@
     )+ };
 }
 
+/// Implements [`Mapping`]
 macro_rules! impl_mapping {
     ($($n:literal)+) => { $(
         impl<'py, const O: u32> Mapping<Loc<$n, f64>> for DolfinxPyFunction_f64<'py, $n, O, 1> {
@@ -342,6 +357,17 @@
 */
 
 impl<'py, const N: u32, const O: u32, const D: u32> DolfinxPyFunction_f64<'py, N, O, D> {
+    fn new_cxx(&mut self, new: *mut Function_f64) {
+        if let None = self.u {
+            if self.cxx != std::ptr::null_mut() {
+                unsafe { ffi::drop_Function_f64(self.cxx) };
+            }
+        }
+        self.cxx = new;
+        self.u = None;
+        self.owned = true;
+    }
+
     /// If owned, call `f_valid` with mutable data.
     /// If not owned, create a new C++ fem.Function, uninitialised, and call `f_invalid` with its
     /// mutable data, as well as the original immutable data.
@@ -352,17 +378,18 @@
         G: Fn(&mut [f64], &[f64]) -> Z,
     {
         if self.owned {
-            return f_valid(unsafe { ffi::data_mut_Function_f64(self.cxx) });
+            let res = f_valid(unsafe { ffi::data_mut_Function_f64(self.cxx) });
+            unsafe { ffi::scatter_fwd_Function_f64(self.cxx) };
+            return res;
         } else {
             // Create a completely new Function
             let old_data = unsafe { ffi::data_mut_Function_f64(self.cxx) };
             let new = unsafe { ffi::similar_Function_f64(self.cxx) };
             let new_data = unsafe { ffi::data_mut_Function_f64(new) };
             let res = f_invalid(new_data, old_data);
+            unsafe { ffi::scatter_fwd_Function_f64(new) };
             // Invalidate Python reference; we're managing the object now
-            self.cxx = new;
-            self.u = None;
-            self.owned = true;
+            self.new_cxx(new);
             return res;
         }
     }
@@ -376,16 +403,17 @@
         F: Fn(&mut [f64]) -> Z,
     {
         if self.owned {
-            return f(unsafe { ffi::data_mut_Function_f64(self.cxx) });
+            let res = f(unsafe { ffi::data_mut_Function_f64(self.cxx) });
+            unsafe { ffi::scatter_fwd_Function_f64(self.cxx) };
+            return res;
         } else {
             // Create a completely new Function
             let new = unsafe { ffi::similar_Function_f64(self.cxx) };
             let new_data = unsafe { ffi::data_mut_Function_f64(new) };
             let res = f(new_data);
+            unsafe { ffi::scatter_fwd_Function_f64(new) };
             // Invalidate Python reference; we're managing the object now
-            self.cxx = new;
-            self.u = None;
-            self.owned = true;
+            self.new_cxx(new);
             return res;
         }
     }
@@ -404,10 +432,11 @@
         let old_data = unsafe { ffi::data_Function_f64(self.cxx) };
         let new_data = unsafe { ffi::data_mut_Function_f64(new.cxx) };
         f(new_data, old_data);
+        unsafe { ffi::scatter_fwd_Function_f64(new.cxx) };
         new
     }
 
-    /// Similar to [`with_mut_data`], but check compatibility with `other` and also pass its
+    /// Similar to [`Self::with_mut_data`], but check compatibility with `other` and also pass its
     /// immutable data slice as the final argument.
     #[inline]
     fn with_mut_data_and_orig_binop<F, G, Z>(
@@ -424,7 +453,7 @@
             ffi::check_compat_Function_f64(self.cxx, other.cxx).unwrap();
             ffi::data_Function_f64(other.cxx)
         };
-        self.with_mut_data_and_orig(
+        let res = self.with_mut_data_and_orig(
             |data_self_valid| {
                 assert_eq!(data_self_valid.len(), data_other.len());
                 f_valid(data_self_valid, data_other)
@@ -434,10 +463,12 @@
                 assert_eq!(data_orig.len(), data_other.len());
                 f_invalid(data_self_invalid, data_orig, data_other)
             },
-        )
+        );
+        unsafe { ffi::scatter_fwd_Function_f64(self.cxx) };
+        return res;
     }
 
-    /// Similar to [`with_new_data_and_orig`], but check compatibility with `other` and also
+    /// Similar to [`Self::with_new_data_and_orig`], but check compatibility with `other` and also
     /// pass its immutable data slice as the final argument.
     #[inline]
     fn with_new_data_and_orig_binop<F>(
@@ -460,6 +491,7 @@
     }
 }
 
+/// Implements unary operations
 macro_rules! impl_unop {
     ($trait:ident, $fn:ident) => {
         impl<'py, const N: u32, const O: u32, const D: u32> $trait
@@ -500,6 +532,7 @@
 
 impl_unop!(Neg, neg);
 
+/// Implements scalar operations
 macro_rules! impl_scalarop {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
         impl<'py, const N: u32, const O: u32, const D: u32> $trait<f64>
@@ -564,6 +597,7 @@
 impl_scalarop!(Mul, mul, MulAssign, mul_assign);
 impl_scalarop!(Div, div, DivAssign, div_assign);
 
+/// Implements consuming binary operations
 macro_rules! impl_binop_consume {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
         impl<'py, const N: u32, const O: u32, const D: u32> $trait<Self>
@@ -635,6 +669,7 @@
     };
 }
 
+/// Implements binary operations on references
 macro_rules! impl_binop_ref {
     ($trait:ident, $fn:ident) => {
         impl<'py, 'b, const N: u32, const O: u32, const D: u32> $trait<Self>
@@ -665,6 +700,7 @@
     };
 }
 
+/// Implements binary operations, both consuming and on references
 macro_rules! impl_binop {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
         impl_binop_consume!($trait, $fn, $trait_assign, $fn_assign);
@@ -677,9 +713,7 @@
 impl_binop!(Add, add, AddAssign, add_assign);
 impl_binop!(Sub, sub, SubAssign, sub_assign);
 
-//use alg_tools::euclidean::Euclidean;
 use alg_tools::linops::{VectorSpace, AXPY};
-//use alg_tools::norms::{Dist, HasDual, Norm, Normed, L2};
 use alg_tools::self_ownable;
 
 self_ownable!(DolfinxPyFunction_f64<'py, N, O, D> where 'py, const N: u32, const O: u32, const D: u32);