Fix UWexpression parameter display and tensor assignment

src/underworld3/function/expressions.py

@@ -106,7 +106,12 @@ def _unwrap_atom(atom, mode='nondimensional'):
                 except Exception:
                     pass
             # Recursively unwrap to get inner expression
-            return atom.sym
+            inner = atom.sym
+            # If inner is a UWQuantity (not SymPy-compatible), resolve it
+            # to a plain value so subs() doesn't fail with SympifyError.
+            if isinstance(inner, UWQuantity) and not isinstance(inner, UWexpression):
+                return _unwrap_atom(inner, mode)
+            return inner
         else:  # symbolic
             return atom.sym
 

src/underworld3/function/unit_conversion.py

@@ -123,11 +123,18 @@ def has_units(obj):
         return True
 
     return False
-# DEPRECATED: Old get_units() from function.unit_conversion has been removed.
-# Use uw.get_units() from the units module instead, which provides the unified,
-# high-level API for extracting units from any object type. The units module
-# version delegates to _extract_units_info() which provides the same smart
-# extraction strategy (prioritizing variables over atoms, avoiding blind tree-walking).
+
+
+def get_units(obj):
+    """Extract unit information from an object.
+
+    Delegates to ``underworld3.units.get_units`` which is the canonical
+    implementation.  This wrapper exists because compiled Cython modules
+    (e.g. ``ckdtree.pyx``) import from this location.
+    """
+    from underworld3.units import get_units as _get_units
+
+    return _get_units(obj)
 
 
 def compute_expression_units(expr):

src/underworld3/maths/tensors.py

@@ -146,6 +146,8 @@ def _unscaled_matrix_to_rank2(
     for I in range(imapping[dim][0]):
         indices = imapping[dim][1]
         i, j = indices[I]
+        v_ij[i, j] = V_I[I]
+        v_ij[j, i] = V_I[I]
 
     return v_ij
 
@@ -167,11 +169,22 @@ def _unscaled_matrix_to_rank4(
         i, j = indices[I]
         for J in range(imapping[dim][0]):
             k, l = indices[J]
+
+            val = C_IJ[I, J]
+
+            # Wrap values that have __getitem__ (e.g. UWexpression from
+            # MathematicalMixin) to prevent SymPy's _setter_iterable_check
+            # from rejecting them during NDimArray assignment.
+            if hasattr(val, "__getitem__") and not isinstance(
+                val, (sympy.MatrixBase, sympy.NDimArray)
+            ):
+                val = sympy.Mul(sympy.S.One, val, evaluate=False)
+
             # C_IJ -> C_ijkl -> C_jilk -> C_ij_lk -> C_jikl (Symmetry)
-            c_ijkl[i, j, k, l] = C_IJ[I, J]
-            c_ijkl[j, i, k, l] = C_IJ[I, J]
-            c_ijkl[i, j, l, k] = C_IJ[I, J]
-            c_ijkl[j, i, l, k] = C_IJ[I, J]
+            c_ijkl[i, j, k, l] = val
+            c_ijkl[j, i, k, l] = val
+            c_ijkl[i, j, l, k] = val
+            c_ijkl[j, i, l, k] = val
 
     return c_ijkl
 

src/underworld3/meshing/surfaces.py

@@ -269,15 +269,13 @@ def _interpolate_to_proxy(self) -> None:
         if self._proxy is None:
             self._create_proxy()
 
-        # Get local mesh node coordinates
-        mesh_coords = self.surface.mesh.X.coords
-        if hasattr(mesh_coords, 'magnitude'):
-            mesh_coords = mesh_coords.magnitude
-        elif hasattr(mesh_coords, '__array__'):
-            mesh_coords = np.asarray(mesh_coords)
-
-        # Get surface vertex data
-        surface_coords = self.surface.vertices
+        # Get local mesh node coordinates in model (internal) space
+        # to match surface vertex coordinates
+        mesh_coords = np.asarray(self.surface.mesh._coords)
+
+        # Get surface vertex data in model (internal) space — bypass the
+        # dimensionalising output gateway so coordinates match mesh._coords.
+        surface_coords = np.array(self.surface._pv_mesh.points) if self.surface._pv_mesh is not None else None
         surface_values = self.data
 
         # For 2D surfaces, use only x,y components for KDTree
@@ -458,24 +456,72 @@ def symbol(self, value: str) -> None:
                 self._distance_var = None
                 self._distance_stale = True
 
+    def _dimensionalise_coords(self, coords: np.ndarray) -> np.ndarray:
+        """Apply dimensional scaling to internal model coordinates.
+
+        Follows the same gateway pattern as ``mesh.X.coords``: internal
+        storage is in model (ND) space; user-facing properties return
+        physical (dimensional) coordinates when the units system is active.
+
+        Returns the original array unchanged when no units are configured.
+        """
+        if coords is None:
+            return None
+
+        cs = getattr(self.mesh, "CoordinateSystem", None)
+        if cs is not None and getattr(cs, "_scaled", False):
+            coords = coords * cs._length_scale
+
+        if self.mesh.units is not None:
+            from underworld3.utilities.unit_aware_array import UnitAwareArray
+
+            return UnitAwareArray(coords, units="meter")
+
+        return coords
+
     @property
     def control_points(self) -> Optional[np.ndarray]:
-        """(N, 3) array of control points defining the surface."""
-        return self._control_points
+        """(N, 3) array of control points defining the surface.
+
+        Returns coordinates in physical (dimensional) units when the
+        units system is active, matching the ``mesh.X.coords`` convention.
+        Internally, points are stored in model (non-dimensional) space.
+        """
+        if self._control_points is None:
+            return None
+        return self._dimensionalise_coords(np.array(self._control_points))
 
     def set_control_points(self, points: np.ndarray) -> None:
         """Set control points and mark discretization as stale.
 
+        Coordinates are stored internally in model (non-dimensional) space,
+        matching the mesh's internal coordinate representation (``mesh._coords``).
+
         Args:
-            points: (N, 2) or (N, 3) array of points. Can include Pint units,
-                   which will be converted using the model's scaling system.
-                   For 2D points, z-coordinate is set to 0.
+            points: (N, 2) or (N, 3) array of points in one of these forms:
+
+                - **Raw numpy array**: Assumed to be in model coordinates
+                  (same space as ``mesh._coords``). When units are active,
+                  this means nondimensional coordinates.
+                - **Pint Quantity / UnitAwareArray**: Automatically converted
+                  to model coordinates via the scaling system.
+
+                For 2D points, a z=0 column is appended automatically.
         """
-        # Handle unit conversion using standard scaling system (ndim)
-        # This works with both legacy (get_coefficients) and modern patterns
-        if hasattr(points, 'magnitude'):
-            # Points have Pint units - convert using scaling system
-            points = uw.scaling.non_dimensionalise(points)
+        # Handle unit conversion: Pint Quantity or UnitAwareArray → model coords
+        # Use the mesh's coordinate system scale factor for consistency with
+        # the output gateway (_dimensionalise_coords).
+        if hasattr(points, "magnitude"):
+            cs = getattr(self.mesh, "CoordinateSystem", None)
+            if cs is not None and getattr(cs, "_scaled", False):
+                # Convert to base SI (meters) then divide by scale factor
+                if hasattr(points, "to_base_units"):
+                    points = np.asarray(points.to_base_units().magnitude) / cs._length_scale
+                else:
+                    points = np.asarray(points.magnitude) / cs._length_scale
+            else:
+                # No scaling active — just strip units
+                points = np.asarray(points.magnitude)
 
         points = np.asarray(points)
         if points.ndim != 2 or points.shape[1] not in (2, 3):
@@ -500,10 +546,15 @@ def _mark_all_proxies_stale(self) -> None:
 
     @property
     def vertices(self) -> Optional[np.ndarray]:
-        """(N, 3) array of vertex positions."""
+        """(N, 3) array of vertex positions.
+
+        Returns coordinates in physical (dimensional) units when the
+        units system is active, matching the ``mesh.X.coords`` convention.
+        Internally, vertices are stored in model (non-dimensional) space.
+        """
         if self._pv_mesh is None:
             return None
-        return np.array(self._pv_mesh.points)
+        return self._dimensionalise_coords(np.array(self._pv_mesh.points))
 
     @property
     def n_vertices(self) -> int:
@@ -783,12 +834,12 @@ def _compute_distance_field(self) -> None:
                 varsymbol=f"d_{{{self._symbol}}}",
             )
 
-        # Get mesh coordinates
-        coords = self.mesh.X.coords
-        if hasattr(coords, 'magnitude'):
-            coords = coords.magnitude
-        elif hasattr(coords, '__array__'):
-            coords = np.asarray(coords)
+        # Get mesh coordinates in model (internal) space.
+        # Must use raw model coordinates (mesh._coords) rather than mesh.X.coords
+        # because mesh.X.coords returns dimensional/scaled coordinates when units
+        # are active, while surface points (self._pv_mesh) are stored in the same
+        # coordinate space used to create them — typically model coordinates.
+        coords = np.asarray(self.mesh._coords)
 
         if self.is_2d:
             # 2D: Use geometry_tools for signed distance to polyline
@@ -1385,12 +1436,9 @@ def compute_distance_field(
                 )
             distance_var = self._distance_var
 
-        # Get mesh coordinates
-        coords = mesh.X.coords
-        if hasattr(coords, 'magnitude'):
-            coords = coords.magnitude
-        elif hasattr(coords, '__array__'):
-            coords = np.asarray(coords)
+        # Get mesh coordinates in model (internal) space
+        # to match surface point coordinates
+        coords = np.asarray(mesh._coords)
 
         pv_mesh = pv.PolyData(coords)
 
@@ -1439,11 +1487,11 @@ def transfer_normals(
         for name, surface in self.surfaces.items():
             surface._ensure_discretized()
 
-        # Get query coordinates
+        # Get query coordinates in model (internal) space
+        # to match surface point coordinates
         if coords is None:
-            coords = mesh.X.coords
-
-        if hasattr(coords, 'magnitude'):
+            coords = np.asarray(mesh._coords)
+        elif hasattr(coords, 'magnitude'):
             coords = coords.magnitude
         elif hasattr(coords, '__array__'):
             coords = np.asarray(coords)

src/underworld3/utilities/_api_tools.py

@@ -290,31 +290,15 @@ def __set__(self, obj, value):
         from ..function.quantities import UWQuantity
         from ..function.expressions import UWexpression
 
-        # Special case: UWexpression assignment (update ._sym AND copy metadata)
-        # UWexpression has ._sym attribute for symbolic value
+        # Special case: UWexpression assignment — store the expression itself
+        # as a symbolic reference (not its inner ._sym value).
+        # This preserves:
+        #   - Symbolic display: parameter shows η₀ not "1e21 Pa.s"
+        #   - Lazy evaluation: if user later changes expression value, parameter updates
+        #   - Unit chain: container.units → value.units → value._sym.units (UWQuantity)
+        # At unwrap/JIT time the chain is followed automatically.
         if isinstance(value, UWexpression):
-            # Update the symbolic value AND copy unit metadata
-            # The expression object (being a SymPy Symbol) preserves identity
-            # while ._sym contains the value for JIT substitution
-            expr.sym = value._sym  # Update substitution value
-
-            # Copy unit metadata from UWexpression
-            if hasattr(value, '_pint_qty'):
-                expr._pint_qty = value._pint_qty
-            if hasattr(value, '_has_pint_qty'):
-                expr._has_pint_qty = value._has_pint_qty
-            if hasattr(value, '_dimensionality'):
-                expr._dimensionality = value._dimensionality
-            if hasattr(value, '_custom_units'):
-                expr._custom_units = value._custom_units
-            if hasattr(value, '_has_custom_units'):
-                expr._has_custom_units = value._has_custom_units
-            if hasattr(value, '_model_registry'):
-                expr._model_registry = value._model_registry
-            if hasattr(value, '_model_instance'):
-                expr._model_instance = value._model_instance
-            if hasattr(value, '_symbolic_with_units'):
-                expr._symbolic_with_units = value._symbolic_with_units
+            expr.sym = value  # Store symbolic reference, not inner value
 
         # Special case: Plain UWQuantity (not UWexpression) - has ._value and ._pint_qty
         # These are simple numbers with units, not symbolic expressions