NVIDIA · peterdsharpe · Sep 8, 2026 · Sep 8, 2026 · Sep 9, 2026 · Sep 9, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -25,12 +25,42 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   `PHYSICSNEMO_DIST_TIMEOUT_S`; unset or empty configuration keeps PyTorch's
   backend default. Invalid timeouts are rejected before initialization state
   changes, allowing corrected configuration to be retried.
+- `MeshToDomainMesh` in `cell_centroids` mode records each source cell's
+  complete effective measure on the interior under the mesh-owned
+  `_effective_measure` point-data key, so
+  integrals and weighted losses over the query points remain possible after
+  the cells are gone.
 - Unified external aero recipe: `NonDimensionalizeByMetadata` gains
   `scale_geometry` so chained instances scale the geometry once; inference
   re-dimensionalizes with the field maps of every instance.
 
 ### Changed
 
+- Mesh integration uses a shared `_effective_measure` field for complete cell
+  and point measures. Cell measures fall back to geometry; point measures are
+  explicit and independent of connectivity. `Mesh.integrate_samples` evaluates
+  point quadrature separately from existing cell and vertex-field integration.
+  Sampling, centroid conversion, geometric transformations, subdivision and
+  GLOBE use the mesh-owned measure API. Point measures carry their represented
+  dimension so geometric scaling preserves their physical units.
+
+  **Migration from 2.2.x:** meshes saved with `cell_data["_measure_weights"]`
+  must be regenerated or converted once before integration:
+
+  ```python
+  from physicsnemo.mesh.calculus import set_cell_measures
+
+  if "_measure_weights" in mesh.cell_data:
+      weights = mesh.cell_data.pop("_measure_weights")
+      set_cell_measures(mesh, mesh.cell_areas * weights)
+  ```
+
+  Replace `compose_measure_weights` calls with `scale_measures`. Consumers
+  should read complete measures with `cell_measures` instead of multiplying
+  `cell_areas` by `cell_measure_weights`. Update stored-field mappings from
+  `cell_data._measure_weights` to `cell_data._effective_measure` and remove any
+  subsequent multiplication by geometric areas.
+
 - `Mesh.slice_points` picks its cell-remapping algorithm by mesh shape: the
   full-mesh lookup table as before, or a binary search over the kept ids when the
   mesh has far more points than cell-vertex entries (a reader keeping a block of
@@ -62,6 +92,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- The unified external aero recipe's force documentation and subsampling warning
+  explain how effective measures compensate for retained-area shrinkage and
+  when sampling or sample-dependent moment origins can introduce bias. Saved
+  inference outputs retain explicit query measures and scale them with geometry
+  so the exported fields remain integrable in physical coordinates.
 - Fixes mesh dtype handling: preserves integer-coordinate precision, normalizes
   connectivity safely, and rejects integer `.to()` casts. Floating/complex casts
   preserve the source mesh.

@@ -50,6 +50,43 @@ Key Operators
 - **Curl**: :math:`\operatorname{curl}(\mathbf{v})` (vector :math:`\to` vector, 3D only)
 - **Laplacian**: :math:`\Delta\varphi` (scalar :math:`\to` scalar, Laplace-Beltrami)
 
+Effective measures and integration
+----------------------------------
+
+The reserved ``_effective_measure`` field stores one complete integration
+measure per cell or point in ``cell_data`` or ``point_data``. Read these with
+``cell_measures(mesh)`` or ``point_measures(mesh)``; they already include any
+geometric contribution and sampling correction.
+
+* ``mesh.integrate(field, data_source="cells")`` integrates piecewise-constant
+  cell values; ``data_source="points"`` integrates piecewise-linear vertex
+  values over the same cells. Both use cell measures, which default to the
+  geometric simplex measures.
+* ``mesh.integrate_samples(field)`` sums independent point samples times their
+  point measures, regardless of connectivity. Explicit point measures are
+  required. For counting measure, use an ordinary sum.
+
+Use ``set_cell_measures`` or ``set_point_measures`` to assign measures. Point
+measures require their represented ``dimension``: 0 for counting, 1 for length,
+2 for area, or 3 for volume. ``scale_measures`` multiplies existing measures by
+a scalar or per-entity factor, such as an inverse sampling probability. Raw
+slicing does not apply a sampling correction.
+
+Converting cells to centroid samples transfers their measures automatically:
+
+.. code:: python
+
+    queries = mesh.to_point_cloud(point_source="cell_centroids")
+    values = queries.points[:, 0]  # Integrate f(x, ...) = x.
+    integral = queries.integrate_samples(values)
+
+For vertex quadrature, ``lumped_point_measures(mesh)`` distributes each cell's
+measure equally among its vertices without modifying the mesh. For finite
+fields, these weights reproduce piecewise-linear integration.
+
+See :doc:`transformations` for supported geometric changes and explicit
+preservation of reference measures.
+
 API Reference
 -------------
 

@@ -36,7 +36,11 @@
 
 from physicsnemo.experimental.utils import CachedPreprocessingDataset
 from physicsnemo.mesh import Mesh
-from physicsnemo.mesh.calculus.measure import compose_measure_weights
+from physicsnemo.mesh.calculus.measure import (
+    cell_measures,
+    scale_measures,
+    set_cell_measures,
+)
 from physicsnemo.mesh.io import from_pyvista
 from physicsnemo.mesh.primitives.planar import structured_grid
 from physicsnemo.mesh.projections import embed
@@ -250,9 +254,8 @@ def subsample_mesh(
           areas that approximate the surface Voronoi diagram of the subsampled
           centroids, at the cost of an O(N) nearest-neighbour pass over the
           full mesh each time a sample is loaded.
-        * **Uniform** (``voronoi=False``): records a single global
-          measure weight (see :mod:`physicsnemo.mesh.calculus.measure`) so
-          that the effective sampled area total matches the full surface.
+        * **Uniform** (``voronoi=False``): rescales the retained effective
+          measures so that their total matches the full surface.
           Much faster, but every subsampled cell gets the same factor
           regardless of local density.
 
@@ -269,7 +272,7 @@ def subsample_mesh(
 
         Returns:
             Mesh with ``n_cells`` cells and a corrected integration measure
-            (uniform: measure weights; Voronoi: area/normal cache).
+            (uniform: reweighted measures; Voronoi: cluster measures and normals).
         """
         if n_cells <= 0:
             raise ValueError(f"{n_cells=!r} must be positive.")
@@ -281,26 +284,22 @@ def subsample_mesh(
         )
 
         if geometry_only:
+            measures = cell_measures(boundary)
             boundary = boundary.with_data(point_data={}, cell_data={}, global_data={})
+            set_cell_measures(boundary, measures)
 
         if voronoi:
-            ### The Voronoi path replaces both areas AND normals with
-            ### locally-accumulated cluster values -- reconstructing the local
-            ### measure and normals, not just correcting the measure -- so it
-            ### keeps the geometry-cache override rather than measure weights.
+            ### Cluster measures are explicit represented areas; the geometric
+            ### cell-area cache continues to describe the retained triangles.
+            ### Cluster normals retain the established boundary-input convention.
             partition = partition_cells(mesh, seeds=boundary.cell_centroids)
-            boundary._cache["cell", "areas"] = partition.cluster_areas
+            set_cell_measures(boundary, partition.cluster_areas)
             boundary._cache["cell", "normals"] = partition.cluster_normals
         else:
-            ### Uniform measure correction via measure weights (see
-            ### physicsnemo.mesh.calculus.measure): one global factor makes the
-            ### effective sampled area total match the full surface exactly
-            ### for this realization (self-normalized / Hajek estimator).
-            ### GLOBE consumes cell_areas * measure_weights, so this is
-            ### numerically identical to the previous area-cache override
-            ### while keeping cell_areas purely geometric.
-            total_area = mesh.cell_areas.sum()
-            compose_measure_weights(boundary, total_area / boundary.cell_areas.sum())
+            ### Normalize the retained measures to preserve the full represented
+            ### area, including any earlier corrections (Hajek estimator).
+            total_area = cell_measures(mesh).sum()
+            scale_measures(boundary, total_area / cell_measures(boundary).sum())
 
         return boundary
 

@@ -38,7 +38,8 @@
 
 The surface integral is evaluated with the mesh quadrature utility
 :meth:`physicsnemo.mesh.Mesh.integrate` (cell-data / P0 rule:
-:math:`\\int_S f\\,dA = \\sum_c f_c\\,A_c`).
+:math:`\\int_S f\\,dA = \\sum_c f_c\\,\\mu_c`, where :math:`\\mu_c` is the
+complete effective measure returned by ``cell_measures(mesh)``).
 
 The coefficient vectors are then projected onto an orthonormal
 (drag, lift, side) triad built from the per-sample freestream direction
@@ -66,16 +67,23 @@
   L_ref`` to integrate on a physical-scale surface; areas and moment
   arms are translation-invariant, so the lost ``CenterMesh`` offset does
   not affect forces (and only shifts the moment reference for moments).
-- **Full surface resolution.** The quadrature covers exactly the cells
-  present on the ``vehicle`` mesh. If the pipeline subsampled the
-  surface (``sampling_resolution`` below the mesh's cell count), the
-  integral covers only the kept cells and every coefficient shrinks by
-  roughly the kept-to-total area fraction -- for predicted and reference
-  values alike, so the *comparison* stays meaningful but the magnitudes
-  do not. ``ForceContext.coefficients``'s 1:1 points/cells contract
-  check cannot detect this (a subsampled surface still satisfies it);
-  ``infer.py`` warns when a vehicle's cell count sits at the
-  ``sampling_resolution`` cap.
+- **Subsampled surfaces.** ``SubsampleMesh`` multiplies each retained cell's
+  effective measure by ``n_before / n_kept``. ``Mesh.integrate`` uses that
+  complete measure to compensate for the retained-area shrinkage. This gives
+  an unbiased Horvitz--Thompson estimate when every cell has inclusion
+  probability ``n_kept / n_before`` and its field value and physical moment
+  reference are fixed independently of the sample.
+  The large-population ``poisson_gap`` sampler is approximate, so this
+  guarantee does not apply to every ``SubsampleMesh`` path. Predictions
+  that depend on the sampled geometry can introduce additional bias.
+  Likewise, centering after reader-level subsampling changes the physical
+  moment origin between samples; reweighting cannot correct that frame change.
+  Subsampled coefficients can have both sampling noise and bias and should be
+  checked for convergence with surface resolution.
+  ``ForceContext.coefficients``'s 1:1 points/cells contract check cannot
+  tell an exact full-surface integral from such an estimate; ``infer.py``
+  warns when a vehicle's cell count sits at the ``sampling_resolution``
+  cap.
 """
 
 from dataclasses import dataclass
@@ -153,7 +161,8 @@ def force_moment_coefficients(
     Args:
         vehicle: Triangulated surface mesh (codimension-1) carrying the
             body. Its cells define the quadrature; ``cell_normals`` and
-            ``cell_areas`` are taken from this mesh.
+            the complete effective cell measures (including sampling corrections)
+            are taken from this mesh.
         pressure_coeff: Per-cell pressure coefficient :math:`C_p`, shape
             ``(n_cells,)`` (a trailing singleton dim, e.g. ``(n_cells, 1)``,
             is flattened internally). Must align 1:1 with ``vehicle`` cells.
@@ -189,7 +198,7 @@ def force_moment_coefficients(
     ### Surface traction coefficient t = -Cp n + Cf (force per area / q_inf).
     traction = -cp.unsqueeze(-1) * normals + cf  # (C, 3)
 
-    ### Quadrature: Mesh.integrate sums field_c * area_c over the surface.
+    ### Quadrature: Mesh.integrate sums field_c * effective_measure_c.
     force = mesh.integrate(traction, data_source="cells")  # (3,)
     arm = centroids - moment_center.to(device=device, dtype=dtype)  # (C, 3)
     moment = mesh.integrate(

@@ -281,13 +281,15 @@ def attach_and_save(
     Writes ``pred_<name>`` and ``true_<name>`` onto a copy of the
     interior's ``point_data`` (the training-space target fields are
     dropped to avoid ambiguity with their physical ``true_<name>``
-    counterparts; non-target inputs like ``sdf`` are kept). The result is
-    saved with :meth:`DomainMesh.save` as a native ``.pdmsh`` tree.
+    counterparts; non-target inputs like ``sdf`` are kept). Explicit point
+    measures are retained and follow geometric rescaling, so the saved sample
+    can still be integrated in physical coordinates. The result is saved with
+    :meth:`DomainMesh.save` as a native ``.pdmsh`` tree.
 
     When *rescale_geometry* is set and ``L_ref`` is available, every mesh
     in the domain is scaled by ``L_ref`` to recover physical-scale
-    coordinates (``Mesh.scale`` leaves ``point_data`` untouched, so the
-    attached fields are not affected).
+    coordinates. ``Mesh.scale`` leaves ordinary ``point_data`` untouched;
+    effective measures scale with their represented dimension.
     """
     if rescale_geometry and "L_ref" in domain.global_data:
         L_ref = domain.global_data["L_ref"]
@@ -567,7 +569,7 @@ def main(cfg: DictConfig) -> None:
     totals: dict[str, float] = {k: 0.0 for k in metric_calculator.expected_keys()}
     count = 0
     sampling_cap = cfg.get("sampling_resolution", None)
-    truncation_warned = False
+    subsampling_warned = False
     for i, idx in enumerate(sampler):
         sample = dataset[idx]
         domain, metadata = sample
@@ -598,23 +600,25 @@ def main(cfg: DictConfig) -> None:
             )
             if sample_forces is not None:
                 force_acc.update(*sample_forces)
-                ### Force magnitudes are only physical at full surface
-                ### resolution (see forces.py): a vehicle cell count
-                ### sitting exactly at the subsample cap means the surface
-                ### was almost certainly truncated by the reader.
+                ### A vehicle cell count sitting exactly at the subsample
+                ### cap means the surface was almost certainly subsampled.
+                ### Effective measures compensate for retained-area shrinkage;
+                ### sampling and moment-frame caveats remain (see forces.py).
                 if (
-                    not truncation_warned
+                    not subsampling_warned
                     and sampling_cap is not None
                     and domain.boundaries["vehicle"].n_cells == sampling_cap
                 ):
                     logger.warning(
                         f"Vehicle surface has exactly sampling_resolution="
                         f"{sampling_cap} cells, so it was likely subsampled; "
-                        f"integrated force/moment coefficients cover only the "
-                        f"kept cells and their magnitudes are not physical. "
-                        f"Raise `sampling_resolution` for absolute CD/CL/CM."
+                        f"integrated force/moment coefficients are estimates "
+                        f"from weighted kept cells and may have sampling "
+                        f"noise or bias, including from a sample-dependent "
+                        f"moment origin. Check convergence by increasing "
+                        f"`sampling_resolution`; see forces.py for assumptions."
                     )
-                    truncation_warned = True
+                    subsampling_warned = True
 
         ### Re-dimensionalize predictions + reference to physical units,
         ### then write them back onto the DomainMesh.

@@ -40,6 +40,7 @@
 from omegaconf import OmegaConf
 
 from physicsnemo.mesh import Mesh
+from physicsnemo.mesh.calculus.measure import scale_measures
 
 
 def _closed_tetrahedron() -> Mesh:
@@ -95,6 +96,34 @@ def test_uniform_shear_gives_drag_equal_to_c_times_area():
     assert abs(res["CL"]) < 1e-4 and abs(res["CS"]) < 1e-4
 
 
+@pytest.mark.parametrize("length_scale", [1.0, 3.0])
+def test_subsampled_surface_with_measure_weights_recovers_full_surface_drag(
+    length_scale,
+):
+    """Measure weights make the integral over kept cells estimate the full one."""
+    full = _closed_tetrahedron()
+    retained = full.slice_cells(torch.tensor([0, 1]))
+    scale_measures(retained, full.n_cells / retained.n_cells)
+
+    c = 2.0
+    cf = torch.tensor([[c, 0.0, 0.0]]).repeat(retained.n_cells, 1)
+    res = forces.force_moment_coefficients(
+        retained,
+        torch.zeros(retained.n_cells),
+        cf,
+        **{**_COMMON, "length_scale": length_scale},
+    )
+
+    ### A regular tetrahedron's faces have equal area, so two of four cells
+    ### weighted by 4/2 reproduce the full-surface value exactly; without
+    ### the weights the integral would be half of it.
+    expected = c * float(full.cell_areas.sum()) * length_scale**2
+    assert res["CD"] == pytest.approx(expected, abs=1e-3)
+    assert c * float(retained.cell_areas.sum()) * length_scale**2 == pytest.approx(
+        expected / 2, abs=1e-3
+    )
+
+
 def test_uniform_shear_moment_about_offset_center():
     """Uniform traction + offset moment center -> analytic yaw moment.
 

@@ -29,13 +29,15 @@
 from pathlib import Path
 
 import infer
+import pytest
 import torch
 from conftest import make_surface_domain_mesh, make_volume_domain_mesh
 from nondim import NonDimensionalizeByMetadata, freestream_scales
 from omegaconf import OmegaConf
 from tensordict import TensorDict
 
 from physicsnemo.mesh import DomainMesh
+from physicsnemo.mesh.calculus.measure import point_measures, set_point_measures
 
 _RECIPE = Path(__file__).resolve().parent.parent
 _DATASETS = _RECIPE / "datasets"
@@ -292,3 +294,29 @@ def test_attach_and_save_rescale_geometry_scales_points(tmp_path):
 
     reloaded = DomainMesh.load(str(out_path))
     assert torch.allclose(reloaded.interior.points, orig_points * l_ref, atol=1e-4)
+
+
+@pytest.mark.parametrize("rescale_geometry", [False, True])
+def test_attach_and_save_preserves_physical_point_measures(tmp_path, rescale_geometry):
+    """Saved measures remain integrable and follow geometry through a round trip."""
+    targets = {"pressure": "scalar", "wss": "vector"}
+    domain = make_surface_domain_mesh(targets, n_cells=16)
+    measures = torch.linspace(0.5, 2.0, domain.interior.n_points)
+    set_point_measures(domain.interior, measures, dimension=2)
+    phys = domain.interior.point_data.select("pressure", "wss")
+    out_path = tmp_path / "m.pdmsh"
+    infer.attach_and_save(
+        domain, phys, phys, targets, out_path, rescale_geometry=rescale_geometry
+    )
+
+    reloaded = DomainMesh.load(str(out_path))
+    factor = domain.global_data["L_ref"] ** 2 if rescale_geometry else 1.0
+    torch.testing.assert_close(point_measures(reloaded.interior), measures * factor)
+    torch.testing.assert_close(
+        reloaded.interior.integrate_samples("pred_pressure"),
+        (phys["pressure"] * measures * factor).sum(),
+    )
+    torch.testing.assert_close(
+        point_measures(reloaded.interior.scale(2.0)), measures * factor * 4
+    )
+    torch.testing.assert_close(point_measures(domain.interior), measures)