From 2194fa9853e3065956a18016fa34574049771dba Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 18 Feb 2026 13:44:48 +1100 Subject: [PATCH 1/3] Fix UWexpression parameter display and tensor assignment in constitutive models MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ExpressionDescriptor.__set__ was extracting the inner value from UWexpression assignments (value._sym), losing the symbolic reference. Parameters displayed as "1e21 Pa.s" instead of the symbolic name (η₀), and lazy evaluation broke. Now stores the UWexpression itself as a symbolic reference, preserving display, lazy evaluation, and the unit delegation chain. Two downstream issues were also fixed: - _unwrap_atom: recursively resolve UWQuantity when found inside UWexpression chain, preventing SympifyError during JIT compilation - _unscaled_matrix_to_rank4: wrap values with __getitem__ (UWexpression from MathematicalMixin) before NDimArray assignment, matching the existing guard in _build_c_tensor. This is the single choke point for all constitutive models after simplify() strips the per-model wrapping. Also fixed pre-existing bug in _unscaled_matrix_to_rank2 where the loop body never assigned values to the output matrix (always returned zeros). Underworld development team with AI support from Claude Code --- src/underworld3/function/expressions.py | 7 +++++- src/underworld3/maths/tensors.py | 21 ++++++++++++---- src/underworld3/utilities/_api_tools.py | 32 +++++++------------------ 3 files changed, 31 insertions(+), 29 deletions(-) diff --git a/src/underworld3/function/expressions.py b/src/underworld3/function/expressions.py index 41cadce2..c0e87eae 100644 --- a/src/underworld3/function/expressions.py +++ b/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 diff --git a/src/underworld3/maths/tensors.py b/src/underworld3/maths/tensors.py index c776889b..68fa4f90 100644 --- a/src/underworld3/maths/tensors.py +++ b/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 diff --git a/src/underworld3/utilities/_api_tools.py b/src/underworld3/utilities/_api_tools.py index a1c590a1..d4649aab 100644 --- a/src/underworld3/utilities/_api_tools.py +++ b/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 From f664320d2ce5fd7db1a9c6f464a1db35e8ec1ee3 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 19 Feb 2026 14:46:00 +1100 Subject: [PATCH 2/3] Restore get_units() in unit_conversion.py for ckdtree compatibility The function was removed during a cleanup but ckdtree.pyx (compiled Cython) still imports it. Added a thin delegate to units.get_units(). Underworld development team with AI support from Claude Code --- src/underworld3/function/unit_conversion.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/underworld3/function/unit_conversion.py b/src/underworld3/function/unit_conversion.py index 4089bc2e..322eb0f6 100644 --- a/src/underworld3/function/unit_conversion.py +++ b/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): From 5da0b6a9f540afc9dae14e7e62eff8e976d4fd3c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 19 Feb 2026 15:00:40 +1100 Subject: [PATCH 3/3] Fix coordinate space handling in Surface and SurfaceVariable Surface stores geometry internally in model (ND) space but user-facing properties (vertices, control_points) must return dimensional coordinates when the units system is active. Changes: - Add _dimensionalise_coords() gateway for Surface output properties - Add set_control_points() input gateway that non-dimensionalises - Fix _interpolate_to_proxy() to use internal ND coordinates for KDTree (was using the now-dimensional .vertices property, causing mismatch with mesh._coords) Underworld development team with AI support from Claude Code --- src/underworld3/meshing/surfaces.py | 122 +++++++++++++++++++--------- 1 file changed, 85 insertions(+), 37 deletions(-) diff --git a/src/underworld3/meshing/surfaces.py b/src/underworld3/meshing/surfaces.py index cd634ea7..e055b47a 100644 --- a/src/underworld3/meshing/surfaces.py +++ b/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)