Materials: a registry and two distributions, and a JIT constant-collision fix

docs/advanced/index.md

@@ -75,6 +75,13 @@ Dynamic remeshing and adaptive refinement strategies.
 
 **[→ Mesh Adaptation](mesh-adaptation.md)**
 
+### Particles: Population Control and Materials
+Keep every cell sampled as the flow deforms the swarm, and declare materials
+on a `MaterialSwarm` so an interface stays where the particles put it — read
+at the integration points, with no level sets in the model script.
+
+**[→ Particle Population and Materials](particle-population-and-materials.md)**
+
 ### Semi-Lagrangian Time Integration (SLCN / SL-BDF2)
 How `AdvDiffusionSLCN` discretizes advection–diffusion in time: the BDF
 time-derivative and Adams-Moulton/θ flux knobs, and how to pair them
@@ -138,6 +145,7 @@ stress-visualisation
 custom-meshes
 curved-boundary-conditions
 mesh-adaptation
+particle-population-and-materials
 semi-lagrangian-time-integration
 eulerian-advection-diffusion
 eulerian-navier-stokes

docs/api/materials.md

@@ -1,27 +1,72 @@
 # Materials
 
+Underworld3 keeps two questions apart: **what** a material is, and **where** it
+is.
+
+*What* is a {py:class}`~underworld3.MaterialRegistry` entry — a name and a
+table of properties, where a property may be a number, a quantity, or a *law*
+such as `eta_0 * sympy.exp(-T.sym[0])`. A registry knows nothing about
+geometry, so it can be built before there is a mesh and shared between models.
+
+*Where* is a **distribution**, and there are two:
+
+- {py:class}`~underworld3.swarm.MaterialSwarm` — carried by particles, so the
+  materials advect with the flow.
+- {py:class}`~underworld3.MaterialRegions` — tied to the mesh, from gmsh
+  physical groups or from a geometric condition. Exact, and needs no particles.
+
+Both present the same face to a solver — `stokes.materials = materials`, which
+sets every constitutive-model parameter the materials declare and the model
+recognises — and both build it from the same partition of unity, one level set
+per material.
+
+See {doc}`../advanced/particle-population-and-materials` for the user guide.
+
 ```{eval-rst}
 .. automodule:: underworld3.materials
    :no-members:
 ```
 
-## Material Property
+## Defining materials
 
 ```{eval-rst}
+.. autoclass:: underworld3.MaterialRegistry
+   :members:
+   :show-inheritance:
+
+.. autoclass:: underworld3.MaterialDefinition
+   :members:
+   :show-inheritance:
+
 .. autoclass:: underworld3.MaterialProperty
    :members:
    :show-inheritance:
 ```
 
-## Material Registry
+## Distributing materials
 
 ```{eval-rst}
-.. autoclass:: underworld3.MaterialRegistry
+.. autoclass:: underworld3.swarm.MaterialSwarm
+   :members:
+   :show-inheritance:
+
+.. autoclass:: underworld3.MaterialRegions
    :members:
    :show-inheritance:
+
+.. autoclass:: underworld3.materials.MaterialDistribution
+   :members:
+   :show-inheritance:
+
+.. autoclass:: underworld3.materials.BoundMaterial
+   :members:
 ```
 
-## Multi-Material Models
+## Different laws per material
+
+A distribution blends per-material *parameter values* into one constitutive
+model. When the materials need genuinely different constitutive **laws** —
+one viscous, one viscoelastic — compose the models instead:
 
 ```{eval-rst}
 .. autoclass:: underworld3.MultiMaterialConstitutiveModel

docs/developer/subsystems/integration-point-variables.md

@@ -69,6 +69,32 @@ symbol onto a `MeshVariable` explicitly with `SNES_Projection`; the
 projection of the field is an ordinary weak form and is exact for data that
 the target space can represent.
 
+
+### The derivative: refused in a weak form, recovered by `evaluate`
+
+An integration-point variable's tabulated gradient is identically zero, so a
+derivative of its symbol would be a silent zero. The two paths are handled
+differently on purpose:
+
+- **Code generation for a weak form** (`utilities/_jitextension.py`,
+  `_no_derivative`) raises. A hidden reconstruction inside a residual would be
+  a per-assembly cost and would decide a discretisation on the user's behalf.
+  The message names the remedy: `proxy_location="cells"`, whose level sets are
+  per-cell polynomials and differentiate directly.
+- **`uw.function.evaluate`** (`function/_function.pyx`,
+  `_integration_point_sources_to_cell_fit`) substitutes any integration-point
+  source appearing under a derivative by a per-cell least-squares fit of its
+  own values, then lets the ordinary derivative machinery run. The fit is
+  allowed to be exactly determined (`nmin = Nb`) because the rule is unisolvent
+  for that degree; the default `Nb + 2` would send every cell to the linear
+  patch and leave the recovered gradient first order.
+
+Measured on `x^2 + 2y` carried at the integration points, the recovered
+gradient converges: 2.4e-3, 6.1e-4, 2.6e-4 at cell sizes 1/5, 1/10, 1/20. The
+direct `"cells"` route (degree 2, fitted from particles) gives 2.4e-7 on the
+same field, because it is exact for a quadratic and nothing is projected
+afterwards.
+
 ## Guards
 
 The field has no gradient (its tabulated derivative is identically zero), so
@@ -335,12 +361,20 @@ nodes and back.
 
 ```python
 swarm = uw.swarm.Swarm(mesh)
-M = uw.swarm.SwarmVariable("M", swarm, 1, proxy_location="integration_points")
+tau = uw.swarm.SwarmVariable("tau", swarm, (2, 2),
+                             proxy_location="integration_points")
 swarm.populate(fill_param=3)
-M.data[:, 0] = ...                                   # per particle
-stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_0 * M.sym[0] + eta_1 * (1 - M.sym[0])
+tau.data[...] = ...                                  # per particle
 ```
 
+For a **material**, this is not the entry point: use `uw.swarm.MaterialSwarm`,
+which owns an `IndexSwarmVariable` whose level sets live at the integration
+points by default and blends the declared properties by the resulting
+partition of unity. See
+{doc}`../../advanced/particle-population-and-materials` for why the direct
+route is not offered — sampling a property field hands the solver an answer
+where it needs a constitutive law.
+
 The reconstruction itself is unchanged (a linear-exact RBF over the nearest
 particles, `rbf_interpolate`); only its target moved. A particle-carried
 material step is reproduced at the integration points with less than half
@@ -350,6 +384,41 @@ the interface (`tests/test_0067_integration_point_proxy.py`).
 has no gradient, so a derivative of the swarm variable's symbol is refused.
 Vector and tensor swarm variables get a multi-component proxy.
 
+### `proxy_sampling`: what each point reads
+
+`proxy_location` is where; `proxy_sampling` is what.
+
+| | `"reconstruct"` (default) | `"share"` |
+|---|---|---|
+| gathers from | the `nnn` nearest particles, by distance | the particles whose nearest integration point *in their own cell* is this one |
+| respects cell walls | no | yes |
+| linear fields | exact | small averaging error |
+| bounded by the particle values | no (overshoots a jump) | yes, it is a mean of them |
+| particles used | the stencil's | all of them, each exactly once |
+
+`"share"` is the cell-restricted Voronoi share
+(`underworld3/utilities/particle_share.py`): `share_assignment` maps each
+particle to one flat index in the cell-major `(ncells, Nq)` layout,
+`share_average` reduces by `np.bincount`. A rule point whose share is empty
+falls back to the nearest particle anywhere on the rank, and the count of
+those is left on `var._share_empty` — persistently non-zero means the swarm is
+too thin for the rule, and `Swarm.repopulate` is the fix.
+
+The assignment needs each particle's owning cell. UW3 swarms are
+`DMSWARM_BASIC`, so PETSc holds no cell id and the locator has to run;
+`Swarm._owning_cells()` caches the result and drops it wherever `_kdtree` is
+dropped, so the share proxies, the population census and anything else
+cell-local pay for one location per step between them. On 32 912 particles /
+242 cells the location is 18.8 ms and the share itself 3.2 ms, against 8.1 ms
+for the `"reconstruct"` path (whose cached operator is geometry-only, so it is
+rebuilt every time the particles move).
+
+There is no `"nearest"` here. Sampling one particle's value whole is the
+material mapping, and materials go through `MaterialSwarm` / its
+`IndexSwarmVariable` (`proxy_sampling="nearest"` there, or `"share"` for
+fractional masks); asking for it on a plain `SwarmVariable` raises and names
+the alternative.
+
 `Lagrangian_Swarm(..., proxy_location="integration_points")` applies the
 same to the fully Lagrangian history: the slots carried on the particles
 are reconstructed at the integration points and the weak form reads them

docs/examples/utilities/intermediate/Ex_Swarm_Material_Index.py

@@ -0,0 +1,149 @@
+# ---
+# jupyter:
+#   jupytext:
+#     formats: py:percent
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#   kernelspec:
+#     display_name: Python 3
+#     language: python
+#     name: python3
+# ---
+
+# %% [markdown]
+"""
+# Materials on a swarm
+
+**PHYSICS:** fluid_mechanics
+**DIFFICULTY:** intermediate
+**PURPOSE:** demonstration
+
+## Description
+
+Two viscosity layers, 1 and 1000, carried by particles and driven from the
+top. The model names its materials and says where they are; it never writes
+a level set, a mask, or a blend.
+
+With the interface on mesh edges the exact velocity is piecewise linear and
+lies in the P2 velocity space, so the only error in the solve is how the
+material is represented. Read at the integration points (the default), where
+each point takes the material of its nearest particle, the problem solves to
+2e-7. Run with `-uw_proxy_location nodes` to watch the nodal level set smear
+the interface across a cell and leave an L2 error of 8e-2.
+"""
+
+# %%
+import numpy as np
+import sympy
+
+import underworld3 as uw
+
+params = uw.Params(
+    uw_proxy_location="integration_points",   # or "nodes", or "cells"
+    uw_proxy_sampling="nearest",              # or "share"
+    uw_cell_size=0.1,
+    uw_eta_top=1000.0,
+    uw_interface=0.5,
+    uw_fill_param=3,
+)
+
+mesh = uw.meshing.UnstructuredSimplexBox(
+    cellSize=params.uw_cell_size, qdegree=2, regular=True
+)
+v = uw.discretisation.MeshVariable("v", mesh, mesh.dim, degree=2)
+p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1)
+
+# %% [markdown]
+"""
+## The materials
+
+`MaterialSwarm` is a `Swarm` that carries materials. Declare them with the
+property names the constitutive model knows, then say where each one is: a
+symbolic condition on the mesh coordinates, a boolean array over the
+particles, or a callable. Every `add` must come before the first read, which
+is what allocates the particles and the level sets.
+"""
+
+# %%
+materials = uw.swarm.MaterialSwarm(
+    mesh,
+    fill_param=params.uw_fill_param,
+    proxy_location=params.uw_proxy_location,
+    proxy_sampling=(
+        params.uw_proxy_sampling
+        if params.uw_proxy_location == "integration_points"
+        else None
+    ),
+)
+materials.add("lower", shear_viscosity_0=1.0, density=3300)
+materials.add("upper", shear_viscosity_0=params.uw_eta_top, density=3400)
+
+materials["upper"] = mesh.X[1] > params.uw_interface
+
+# %% [markdown]
+"""
+## The solve
+
+`stokes.materials = materials` sets every parameter the constitutive model
+recognises — here `shear_viscosity_0`. `density` is not a viscous-model
+parameter, so it stays a blended symbol for the model script to use; this
+problem is driven by the boundary, so it goes unused.
+"""
+
+# %%
+stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+stokes.materials = materials
+
+stokes.add_dirichlet_bc((1.0, 0.0), "Top")
+stokes.add_dirichlet_bc((0.0, 0.0), "Bottom")
+stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left")
+stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right")
+stokes.tolerance = 1e-8
+stokes.solve()
+
+# %% [markdown]
+"""
+## Against the exact layered Couette profile
+
+The error is a global integral rather than a nodal norm, so the number is the
+same however the mesh is partitioned.
+"""
+
+# %%
+h, eta_top = params.uw_interface, params.uw_eta_top
+gradient = 1.0 / (h + (1.0 - h) / eta_top)
+y = mesh.X[1]
+exact = sympy.Piecewise(
+    (gradient * y, y < h),
+    (gradient * h + gradient / eta_top * (y - h), True),
+)
+
+error = uw.maths.Integral(mesh, (v.sym[0] - exact) ** 2).evaluate() ** 0.5
+viscosity = uw.maths.Integral(mesh, materials.shear_viscosity_0).evaluate()
+exact_viscosity = 1.0 * h + eta_top * (1.0 - h)
+
+uw.pprint(
+    f"proxy_location={params.uw_proxy_location} "
+    f"sampling={params.uw_proxy_sampling} fill={params.uw_fill_param}: "
+    f"assembled int(eta) {viscosity:.4f} (exact {exact_viscosity:.4f}) | "
+    f"velocity L2 {error:.3e}"
+)
+
+# %% [markdown]
+"""
+## What the weak form sees
+
+The upper-material mask along a line crossing the interface. At the
+integration points it is a step in the right place; the nodal level set ramps
+across a whole cell.
+"""
+
+# %%
+line = np.column_stack(
+    [np.full(201, 0.5), np.linspace(max(0.0, h - 0.25), min(1.0, h + 0.25), 201)]
+)
+upper = np.asarray(uw.function.evaluate(materials["upper"].mask, line)).reshape(-1)
+uw.pprint(f"mask along the line: {upper.min():+.4f} .. {upper.max():+.4f}")