From 1c87ed1767e450421ec4393147c480c57893d144 Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Tue, 8 Sep 2026 08:49:40 -0400 Subject: [PATCH 1/6] Carry the quadrature measure onto cell-centroid query points MeshToDomainMesh in cell_centroids mode turns each source cell into an interior query point and discards the cells. Any integral over those points (forces, area-weighted losses or metrics) then has no measure to weight by, and after SubsampleMesh the retained cells' measure weights are lost as well. Record cell_measures(mesh) (area times composed measure weights) on the interior under the reserved point_data key TARGET_QUADRATURE_MEASURE_KEY, aligned one-for-one with the centroids. Reject the key as a user target or as a pre-existing input field so it cannot be silently shadowed. Also validate the shape of the reserved per-cell measure-weights field and of tensor factors passed to compose_measure_weights: TensorDict only checks the leading dimension, so a (n_cells, 1) tensor was storable but broadcast wrongly against cell_areas. --- CHANGELOG.md | 5 ++ .../datapipes/transforms/mesh/__init__.py | 2 + .../datapipes/transforms/mesh/transforms.py | 38 +++++++++-- physicsnemo/mesh/calculus/measure.py | 29 +++++++- .../transforms/test_mesh_to_domain_mesh.py | 68 +++++++++++++++++-- .../transforms/test_nested_fields.py | 6 +- test/mesh/calculus/test_measure.py | 18 ++++- 7 files changed, 154 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b252c40d93..7b9d49f472 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 in `physicsnemo.datapipes` through it, so a `"."` in a YAML field name (`"solution.pressure"`) addresses a leaf inside a nested `TensorDict`. Nested `Mesh` data no longer needs to be flattened before use. +- `MeshToDomainMesh` in `cell_centroids` mode records each source cell's + effective measure (area times any composed measure weights) on the interior + under the reserved `point_data` key `TARGET_QUADRATURE_MEASURE_KEY`, so + integrals and weighted losses over the query points remain possible after + the cells are gone. ### Changed diff --git a/physicsnemo/datapipes/transforms/mesh/__init__.py b/physicsnemo/datapipes/transforms/mesh/__init__.py index 0fd65ca02a..b032bb46c9 100644 --- a/physicsnemo/datapipes/transforms/mesh/__init__.py +++ b/physicsnemo/datapipes/transforms/mesh/__init__.py @@ -29,6 +29,7 @@ ) from physicsnemo.datapipes.transforms.mesh.base import MeshTransform from physicsnemo.datapipes.transforms.mesh.transforms import ( + TARGET_QUADRATURE_MEASURE_KEY, CenterMesh, ComputeCellCentroids, ComputeSurfaceNormals, @@ -59,6 +60,7 @@ "SetGlobalField", "NormalizeMeshFields", "MeshToDomainMesh", + "TARGET_QUADRATURE_MEASURE_KEY", "MeshToTensorDict", "RestructureTensorDict", "RandomScaleMesh", diff --git a/physicsnemo/datapipes/transforms/mesh/transforms.py b/physicsnemo/datapipes/transforms/mesh/transforms.py index 73bdb8bb33..63d120d7ac 100644 --- a/physicsnemo/datapipes/transforms/mesh/transforms.py +++ b/physicsnemo/datapipes/transforms/mesh/transforms.py @@ -45,9 +45,16 @@ Mesh, MeshFieldAssociation, ) -from physicsnemo.mesh.calculus.measure import compose_measure_weights +from physicsnemo.mesh.calculus.measure import cell_measures, compose_measure_weights from physicsnemo.nn.functional import weighted_multinomial +### Reserved ``point_data`` key carrying the effective quadrature measure of +### the cell-centroid query points created by :class:`MeshToDomainMesh`. +### Distinct from ``MEASURE_WEIGHTS_KEY``: that is a dimensionless factor on +### source cells, while this is the full geometric measure (area * weight), +### aligned one-for-one with the interior points. +TARGET_QUADRATURE_MEASURE_KEY: str = "_target_quadrature_measure" + @register() class ScaleMesh(MeshTransform): @@ -1003,8 +1010,11 @@ class MeshToDomainMesh(MeshTransform): Names of cell-centered fields on the input mesh to use as prediction targets. They are moved out of the boundary's ``cell_data`` and into ``interior.point_data``. Use with ``interior_points='cell_centroids'``. - If ``None`` (and ``point_data_targets`` is also ``None``), no targets - are placed on the interior. + If ``None`` (and ``point_data_targets`` is also ``None``), no user + targets are placed on the interior. Centroid mode still records each + source cell's effective measure under + :data:`TARGET_QUADRATURE_MEASURE_KEY`, so integrals over the query + points remain possible after the cells are gone. point_data_targets : list[str] or None, default ``None`` Names of vertex-centered fields on the input mesh to use as prediction targets. They are moved out of the boundary's ``point_data`` and into @@ -1086,10 +1096,27 @@ def __init__( ### ``select`` / ``exclude`` below accept the parsed tuple keys. self._cell_data_targets: list[NestedKey] = as_nested_keys(cell_data_targets) self._point_data_targets: list[NestedKey] = as_nested_keys(point_data_targets) + if TARGET_QUADRATURE_MEASURE_KEY in ( + self._cell_data_targets + self._point_data_targets + ): + raise ValueError( + f"{TARGET_QUADRATURE_MEASURE_KEY!r} is reserved for the query " + "measure and cannot be configured as a target." + ) self._interior_points = interior_points self._boundary_name = boundary_name def __call__(self, mesh: Mesh) -> DomainMesh: # type: ignore[override] + for association, data in ( + ("point_data", mesh.point_data), + ("cell_data", mesh.cell_data), + ): + if TARGET_QUADRATURE_MEASURE_KEY in data: + raise ValueError( + f"Input mesh {association} already contains reserved key " + f"{TARGET_QUADRATURE_MEASURE_KEY!r}; rename the field " + "before MeshToDomainMesh." + ) ### v1 supports two diagonal corners: ### (cell_data_targets, interior_points='cell_centroids') ### (point_data_targets, interior_points='vertices') @@ -1123,13 +1150,16 @@ def __call__(self, mesh: Mesh) -> DomainMesh: # type: ignore[override] def _call_cell_centroids(self, mesh: Mesh) -> DomainMesh: ### Build the interior as a point cloud at cell centroids, with target - ### cell_data fields moved into interior.point_data. + ### cell_data fields moved into interior.point_data. The source cells + ### do not exist on the interior, so record their effective measure + ### (area * any composed measure weights) beside the centroids now. require_keys(mesh.cell_data, self._cell_data_targets, what="Target field") interior_point_data = ( mesh.cell_data.select(*self._cell_data_targets) if self._cell_data_targets else TensorDict({}, batch_size=[mesh.n_cells]) ) + interior_point_data[TARGET_QUADRATURE_MEASURE_KEY] = cell_measures(mesh) interior = Mesh( points=mesh.cell_centroids, point_data=interior_point_data, diff --git a/physicsnemo/mesh/calculus/measure.py b/physicsnemo/mesh/calculus/measure.py index 189152dbaf..59bd6c5f8a 100644 --- a/physicsnemo/mesh/calculus/measure.py +++ b/physicsnemo/mesh/calculus/measure.py @@ -80,6 +80,24 @@ MEASURE_WEIGHTS_KEY: str = "_measure_weights" +def _validate_measure_weight_shape( + mesh: "Mesh", weights: torch.Tensor +) -> Float[torch.Tensor, " n_cells"]: + """Enforce the reserved field's one-scalar-per-cell contract. + + ``cell_data``'s batch dimension only checks the leading axis, so a + ``(n_cells, 1)`` tensor is storable but would broadcast wrongly against + ``cell_areas``. + """ + expected = (mesh.n_cells,) + if weights.shape != expected: + raise ValueError( + f"{MEASURE_WEIGHTS_KEY!r} must have shape {expected} " + f"(one scalar per cell), got {tuple(weights.shape)}" + ) + return weights + + def cell_measure_weights(mesh: "Mesh") -> Float[torch.Tensor, " n_cells"]: r"""Per-cell measure weights of *mesh* (ones when none are recorded). @@ -94,7 +112,7 @@ def cell_measure_weights(mesh: "Mesh") -> Float[torch.Tensor, " n_cells"]: return torch.ones( mesh.n_cells, dtype=mesh.points.dtype, device=mesh.points.device ) - return weights + return _validate_measure_weight_shape(mesh, weights) def cell_measures(mesh: "Mesh") -> Float[torch.Tensor, " n_cells"]: @@ -114,7 +132,7 @@ def cell_measures(mesh: "Mesh") -> Float[torch.Tensor, " n_cells"]: weights = mesh.cell_data.get(MEASURE_WEIGHTS_KEY, None) if weights is None: return cell_areas - return cell_areas * weights + return cell_areas * _validate_measure_weight_shape(mesh, weights) def compose_measure_weights(mesh: "Mesh", factor: float | torch.Tensor) -> None: @@ -135,4 +153,11 @@ def compose_measure_weights(mesh: "Mesh", factor: float | torch.Tensor) -> None: This producer's weight contribution; a scalar or a per-cell tensor broadcast against the existing weights. """ + if isinstance(factor, torch.Tensor) and factor.ndim != 0: + expected = (mesh.n_cells,) + if factor.shape != expected: + raise ValueError( + f"measure-weight factor must be a scalar or have shape " + f"{expected}, got {tuple(factor.shape)}" + ) mesh.cell_data[MEASURE_WEIGHTS_KEY] = cell_measure_weights(mesh) * factor diff --git a/test/datapipes/transforms/test_mesh_to_domain_mesh.py b/test/datapipes/transforms/test_mesh_to_domain_mesh.py index 8bcc4765f9..147f85cd20 100644 --- a/test/datapipes/transforms/test_mesh_to_domain_mesh.py +++ b/test/datapipes/transforms/test_mesh_to_domain_mesh.py @@ -19,8 +19,12 @@ import pytest import torch -from physicsnemo.datapipes.transforms.mesh import MeshToDomainMesh +from physicsnemo.datapipes.transforms.mesh import ( + TARGET_QUADRATURE_MEASURE_KEY, + MeshToDomainMesh, +) from physicsnemo.mesh import DomainMesh, Mesh +from physicsnemo.mesh.calculus import MEASURE_WEIGHTS_KEY, cell_measures # --------------------------------------------------------------------------- # Helpers @@ -76,6 +80,27 @@ def _point_cloud_mesh_3d(n_points: int = 5) -> Mesh: ) +def _unequal_triangle_mesh_3d() -> Mesh: + """Two disconnected triangles with areas 0.5 and 2.0 and measure weights.""" + return Mesh( + points=torch.tensor( + [ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [3.0, 0.0, 0.0], + [5.0, 0.0, 0.0], + [3.0, 2.0, 0.0], + ] + ), + cells=torch.tensor([[0, 1, 2], [3, 4, 5]]), + cell_data={ + "C_p": torch.tensor([10.0, 20.0]), + MEASURE_WEIGHTS_KEY: torch.tensor([2.0, 0.25]), + }, + ) + + def _domain_mesh_3d() -> DomainMesh: """A minimal 3-D DomainMesh (interior point cloud + one boundary).""" interior = Mesh( @@ -160,10 +185,32 @@ def test_targets_moved_to_interior_point_data(self): transform = MeshToDomainMesh(cell_data_targets=["C_p", "C_f"]) domain = transform(mesh) interior_keys = set(domain.interior.point_data.keys()) - assert interior_keys == {"C_p", "C_f"} + assert interior_keys == {"C_p", "C_f", TARGET_QUADRATURE_MEASURE_KEY} assert torch.allclose(domain.interior.point_data["C_p"], mesh.cell_data["C_p"]) assert torch.allclose(domain.interior.point_data["C_f"], mesh.cell_data["C_f"]) + def test_target_measure_is_effective_cell_measure_in_cell_order(self): + mesh = _unequal_triangle_mesh_3d() + domain = MeshToDomainMesh(cell_data_targets=["C_p"])(mesh) + boundary = domain.boundaries["vehicle"] + + assert torch.equal( + domain.interior.point_data[TARGET_QUADRATURE_MEASURE_KEY], + cell_measures(mesh), + ) + assert torch.equal( + domain.interior.point_data["C_p"], torch.tensor([10.0, 20.0]) + ) + ### Query i, target i, and measure i all refer to boundary cell i. + assert torch.equal(domain.interior.points, boundary.cell_centroids) + assert torch.equal( + domain.interior.point_data[TARGET_QUADRATURE_MEASURE_KEY], + cell_measures(boundary), + ) + ### The measure lives only on the interior, never on the boundary. + assert TARGET_QUADRATURE_MEASURE_KEY not in boundary.point_data + assert TARGET_QUADRATURE_MEASURE_KEY not in boundary.cell_data + def test_non_target_cell_data_stays_on_boundary(self): mesh = _two_triangle_mesh_3d() transform = MeshToDomainMesh(cell_data_targets=["C_p", "C_f"]) @@ -197,15 +244,28 @@ def test_global_data_passed_through(self): assert "U_inf" in domain.global_data.keys() assert torch.allclose(domain.global_data["U_inf"], mesh.global_data["U_inf"]) - def test_no_targets_yields_empty_interior_point_data(self): + def test_no_targets_yields_only_target_measure(self): mesh = _two_triangle_mesh_3d() transform = MeshToDomainMesh(cell_data_targets=None) domain = transform(mesh) - assert len(domain.interior.point_data.keys()) == 0 + assert set(domain.interior.point_data.keys()) == {TARGET_QUADRATURE_MEASURE_KEY} ### All original cell_data should still be on the boundary. boundary_keys = set(domain.boundaries["vehicle"].cell_data.keys()) assert boundary_keys == {"C_p", "C_f", "normals"} + def test_reserved_measure_key_cannot_be_a_target(self): + with pytest.raises(ValueError, match="reserved"): + MeshToDomainMesh(cell_data_targets=[TARGET_QUADRATURE_MEASURE_KEY]) + + @pytest.mark.parametrize("association", ["point_data", "cell_data"]) + def test_preexisting_reserved_measure_key_is_rejected(self, association): + mesh = _two_triangle_mesh_3d() + getattr(mesh, association)[TARGET_QUADRATURE_MEASURE_KEY] = torch.ones( + mesh.n_points if association == "point_data" else mesh.n_cells + ) + with pytest.raises(ValueError, match="already contains reserved key"): + MeshToDomainMesh(cell_data_targets=["C_p"])(mesh) + def test_custom_boundary_name(self): mesh = _two_triangle_mesh_3d() transform = MeshToDomainMesh(cell_data_targets=["C_p"], boundary_name="airfoil") diff --git a/test/datapipes/transforms/test_nested_fields.py b/test/datapipes/transforms/test_nested_fields.py index bb1adff691..ef63bce357 100644 --- a/test/datapipes/transforms/test_nested_fields.py +++ b/test/datapipes/transforms/test_nested_fields.py @@ -29,6 +29,7 @@ import physicsnemo.datapipes as dp from physicsnemo.datapipes.transforms.mesh import ( + TARGET_QUADRATURE_MEASURE_KEY, ComputeSurfaceNormals, DropMeshFields, MeshToDomainMesh, @@ -253,7 +254,10 @@ def test_missing_or_leaf_prefix_target_raises_key_error(self): def test_nested_target_moved_to_interior(self): mesh = _surface_mesh() domain = MeshToDomainMesh(cell_data_targets=["solution.pMeanTrim"])(mesh) - assert _leaves(domain.interior.point_data) == {("solution", "pMeanTrim")} + assert _leaves(domain.interior.point_data) == { + ("solution", "pMeanTrim"), + TARGET_QUADRATURE_MEASURE_KEY, + } boundary = domain.boundaries["vehicle"] assert ("solution", "pMeanTrim") not in boundary.cell_data assert ("solution", "wssMeanTrim") in boundary.cell_data diff --git a/test/mesh/calculus/test_measure.py b/test/mesh/calculus/test_measure.py index 5ea38f483e..d6b0d228bc 100644 --- a/test/mesh/calculus/test_measure.py +++ b/test/mesh/calculus/test_measure.py @@ -87,11 +87,27 @@ def test_compose_roundtrip_via_reserved_key(self): ) def test_storage_rejects_wrong_shape(self): - ### cell_data's batch dimension enforces the (n_cells,) shape. + ### cell_data's batch dimension rejects a wrong leading dimension. mesh = two_triangles_2d.load() with pytest.raises(RuntimeError): mesh.cell_data[MEASURE_WEIGHTS_KEY] = torch.ones(mesh.n_cells + 1) + def test_reserved_field_rejects_trailing_singleton_dimension(self): + ### TensorDict accepts vector-valued cell data, so the reserved + ### scalar field must enforce its own exact shape. + mesh = two_triangles_2d.load() + mesh.cell_data[MEASURE_WEIGHTS_KEY] = torch.ones(mesh.n_cells, 1) + + with pytest.raises(ValueError, match="one scalar per cell"): + cell_measure_weights(mesh) + with pytest.raises(ValueError, match="one scalar per cell"): + cell_measures(mesh) + + def test_compose_rejects_non_scalar_broadcast_shape(self): + mesh = two_triangles_2d.load() + with pytest.raises(ValueError, match="scalar or have shape"): + compose_measure_weights(mesh, torch.ones(mesh.n_cells, 1)) + def test_weights_survive_slice_cells(self): mesh = make_triangle_strip(6) compose_measure_weights(mesh, torch.arange(1.0, 7.0)) From e39cd39a6e848c9653677c917fe9f980c5f83f05 Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Tue, 8 Sep 2026 08:52:43 -0400 Subject: [PATCH 2/6] Aero recipe: describe subsampled-surface forces as unbiased estimates, drop query measure from saved outputs forces.py and infer.py said that force/moment coefficients from a subsampled vehicle surface "cover only the kept cells" and shrink by the kept-area fraction. That is no longer true: SubsampleMesh records each kept cell's inverse inclusion probability as a measure weight and Mesh.integrate multiplies by it, so the coefficients are unbiased (if noisy) estimates of the full-surface integrals. Correct the module docstring and the once-per-run warning, and pin the behaviour with a test on a subsampled closed surface. MeshToDomainMesh now records the query measure on the interior. It is in training-geometry units and would be stale after rescale_geometry, so attach_and_save drops it alongside the training-space targets. --- CHANGELOG.md | 6 ++++ .../src/forces.py | 27 +++++++------- .../unified_external_aero_recipe/src/infer.py | 35 +++++++++++-------- .../tests/test_forces.py | 21 +++++++++++ .../tests/test_infer.py | 16 +++++++++ 5 files changed, 79 insertions(+), 26 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7b9d49f472..d8eda02c57 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- The unified external aero recipe no longer writes `MeshToDomainMesh`'s + reserved query-measure key into saved inference outputs, and its force + documentation and subsampling warning now describe coefficients from a + subsampled surface as unbiased estimates (the measure weights recorded by + `SubsampleMesh` already enter `Mesh.integrate`) rather than as values shrunk + by the kept-area fraction. - Datapipe transforms, collators, readers, and the unified external aero recipe no longer silently skip or mis-handle nested `TensorDict` fields (membership was tested against top-level `td.keys()`, and diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/forces.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/forces.py index 648953a0f9..9ba8db67d8 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/forces.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/forces.py @@ -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\\,A_c\\,w_c`, where :math:`w_c` is the +cell's measure weight -- one unless the surface was subsampled). The coefficient vectors are then projected onto an orthonormal (drag, lift, side) triad built from the per-sample freestream direction @@ -66,16 +67,17 @@ 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`` records each kept cell's + inverse inclusion probability (``n_before / n_kept``) as its measure + weight, and ``Mesh.integrate`` multiplies by it. Coefficients computed + from a subsampled ``vehicle`` mesh are therefore unbiased estimates of + the full-surface integrals (a Horvitz--Thompson estimator), not values + shrunk by the kept-area fraction -- but they are stochastic, with + sampling noise that grows as fewer cells are kept. + ``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 +155,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 effective cell measures (areas times any measure weights) + 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. diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/infer.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/infer.py index 2ffc05451b..3611c70e99 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/infer.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/infer.py @@ -99,6 +99,7 @@ from physicsnemo import datapipes # noqa: F401 - registers ${dp:...} resolver from physicsnemo.datapipes.keys import as_nested_key, with_leaf_name +from physicsnemo.datapipes.transforms.mesh import TARGET_QUADRATURE_MEASURE_KEY from physicsnemo.distributed import DistributedManager, fused_all_reduce from physicsnemo.mesh import DomainMesh from physicsnemo.utils import load_checkpoint @@ -277,8 +278,9 @@ def attach_and_save( Writes ``pred_`` and ``true_`` onto a copy of the interior's ``point_data`` (the training-space target fields are dropped to avoid ambiguity with their physical ``true_`` - 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; the query + measure ``MeshToDomainMesh`` records is dropped). 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 @@ -295,9 +297,13 @@ def attach_and_save( ### Names may spell nested leaves ("solution.p"); ``key in td`` and ### ``exclude`` resolve them, and the pred_/true_ prefix goes on the ### leaf so the nesting is preserved: ("solution", "pred_p"). + ### The query measure is in training-geometry units and would be stale + ### once the geometry is rescaled, so it is not written either. target_keys = [as_nested_key(n) for n in target_config] - present_targets = [k for k in target_keys if k in interior.point_data] - new_pd = interior.point_data.exclude(*present_targets).clone() + drop_keys = [k for k in target_keys if k in interior.point_data] + if TARGET_QUADRATURE_MEASURE_KEY in interior.point_data: + drop_keys.append(TARGET_QUADRATURE_MEASURE_KEY) + new_pd = interior.point_data.exclude(*drop_keys).clone() for key, val in pred_phys.items(include_nested=True, leaves_only=True): new_pd[with_leaf_name(key, lambda n: f"pred_{n}")] = val for key, val in true_phys.items(include_nested=True, leaves_only=True): @@ -563,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 @@ -594,23 +600,24 @@ 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, + ### so the coefficients are unbiased but noisy estimates of + ### the full-surface integrals (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 the kept cells (unbiased, but with sampling " + f"noise), not exact full-surface integrals. Raise " + f"`sampling_resolution` for exact CD/CL/CM." ) - truncation_warned = True + subsampling_warned = True ### Re-dimensionalize predictions + reference to physical units, ### then write them back onto the DomainMesh. diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_forces.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_forces.py index db8cd5c0a5..83d8e09b90 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_forces.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_forces.py @@ -40,6 +40,7 @@ from omegaconf import OmegaConf from physicsnemo.mesh import Mesh +from physicsnemo.mesh.calculus.measure import compose_measure_weights def _closed_tetrahedron() -> Mesh: @@ -95,6 +96,26 @@ def test_uniform_shear_gives_drag_equal_to_c_times_area(): assert abs(res["CL"]) < 1e-4 and abs(res["CS"]) < 1e-4 +def test_subsampled_surface_with_measure_weights_recovers_full_surface_drag(): + """Measure weights make the integral over kept cells estimate the full one.""" + full = _closed_tetrahedron() + retained = full.slice_cells(torch.tensor([0, 1])) + compose_measure_weights(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 + ) + + ### 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()) + assert res["CD"] == pytest.approx(expected, abs=1e-3) + assert c * float(retained.cell_areas.sum()) == pytest.approx(expected / 2, abs=1e-3) + + def test_uniform_shear_moment_about_offset_center(): """Uniform traction + offset moment center -> analytic yaw moment. diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_infer.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_infer.py index eea6beb52d..f4a6ffca2e 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_infer.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/tests/test_infer.py @@ -35,6 +35,7 @@ from omegaconf import OmegaConf from tensordict import TensorDict +from physicsnemo.datapipes.transforms.mesh import TARGET_QUADRATURE_MEASURE_KEY from physicsnemo.mesh import DomainMesh _RECIPE = Path(__file__).resolve().parent.parent @@ -257,3 +258,18 @@ 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) + + +def test_attach_and_save_drops_query_measure(tmp_path): + """The training-geometry query measure is not written to the artifact.""" + targets = {"pressure": "scalar", "wss": "vector"} + domain = make_surface_domain_mesh(targets, n_cells=16) + domain.interior.point_data[TARGET_QUADRATURE_MEASURE_KEY] = torch.ones( + domain.interior.n_points + ) + 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=True) + + reloaded = DomainMesh.load(str(out_path)) + assert TARGET_QUADRATURE_MEASURE_KEY not in reloaded.interior.point_data From a71c242a86c9e5779dc0c492a8cf05b24947b3ff Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Wed, 9 Sep 2026 09:34:59 -0400 Subject: [PATCH 3/6] Clarify sampling and reference-frame limits of force estimates Signed-off-by: Peter Sharpe --- CHANGELOG.md | 8 ++++---- .../src/forces.py | 20 ++++++++++++------- .../unified_external_aero_recipe/src/infer.py | 13 ++++++------ 3 files changed, 24 insertions(+), 17 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c8c6ac225d..bd6c6d92ba 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -36,10 +36,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - The unified external aero recipe no longer writes `MeshToDomainMesh`'s reserved query-measure key into saved inference outputs, and its force - documentation and subsampling warning now describe coefficients from a - subsampled surface as unbiased estimates (the measure weights recorded by - `SubsampleMesh` already enter `Mesh.integrate`) rather than as values shrunk - by the kept-area fraction. + documentation and subsampling warning now explain how measure weights + compensate for retained-area shrinkage. Exact unbiasedness requires the + correct inclusion probabilities, fixed fields, and a fixed physical moment + origin; approximate sampling and sample-dependent frames can introduce bias. - Datapipe transforms, collators, readers, and the unified external aero recipe no longer silently skip or mis-handle nested `TensorDict` fields (membership was tested against top-level `td.keys()`, and diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/forces.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/forces.py index 9ba8db67d8..e25d6d42ae 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/forces.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/forces.py @@ -67,13 +67,19 @@ 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). -- **Subsampled surfaces.** ``SubsampleMesh`` records each kept cell's - inverse inclusion probability (``n_before / n_kept``) as its measure - weight, and ``Mesh.integrate`` multiplies by it. Coefficients computed - from a subsampled ``vehicle`` mesh are therefore unbiased estimates of - the full-surface integrals (a Horvitz--Thompson estimator), not values - shrunk by the kept-area fraction -- but they are stochastic, with - sampling noise that grows as fewer cells are kept. +- **Subsampled surfaces.** ``SubsampleMesh`` records the cell-count + correction ``n_before / n_kept`` as a measure weight, and + ``Mesh.integrate`` multiplies by it 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; measure weights 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`` diff --git a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/infer.py b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/infer.py index 3611c70e99..7a0da66dc8 100644 --- a/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/infer.py +++ b/examples/cfd/external_aerodynamics/unified_external_aero_recipe/src/infer.py @@ -601,9 +601,9 @@ def main(cfg: DictConfig) -> None: if sample_forces is not None: force_acc.update(*sample_forces) ### A vehicle cell count sitting exactly at the subsample - ### cap means the surface was almost certainly subsampled, - ### so the coefficients are unbiased but noisy estimates of - ### the full-surface integrals (see forces.py). + ### cap means the surface was almost certainly subsampled. + ### Measure weights compensate for retained-area shrinkage; + ### sampling and moment-frame caveats remain (see forces.py). if ( not subsampling_warned and sampling_cap is not None @@ -613,9 +613,10 @@ def main(cfg: DictConfig) -> None: f"Vehicle surface has exactly sampling_resolution=" f"{sampling_cap} cells, so it was likely subsampled; " f"integrated force/moment coefficients are estimates " - f"from the kept cells (unbiased, but with sampling " - f"noise), not exact full-surface integrals. Raise " - f"`sampling_resolution` for exact CD/CL/CM." + 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." ) subsampling_warned = True From 646aafcbb7a9c381a8b8091b207fcb51fb64f005 Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Sat, 19 Sep 2026 16:17:13 -0400 Subject: [PATCH 4/6] Unify complete integration measures across cells and point samples --- CHANGELOG.md | 12 +- docs/api/mesh/calculus.rst | 118 +++++++ .../globe/drivaer/dataset.py | 37 +- physicsnemo/datapipes/readers/mesh.py | 31 +- .../datapipes/transforms/mesh/__init__.py | 2 - .../datapipes/transforms/mesh/transforms.py | 69 ++-- physicsnemo/mesh/calculus/__init__.py | 11 +- physicsnemo/mesh/calculus/integration.py | 24 +- physicsnemo/mesh/calculus/measure.py | 318 ++++++++++------- physicsnemo/mesh/mesh.py | 90 ++++- physicsnemo/mesh/remeshing/_partition.py | 6 +- physicsnemo/mesh/remeshing/_remeshing.py | 10 + physicsnemo/mesh/subdivision/_data.py | 7 + physicsnemo/mesh/subdivision/butterfly.py | 7 +- physicsnemo/mesh/subdivision/linear.py | 7 +- physicsnemo/mesh/subdivision/loop.py | 7 +- physicsnemo/mesh/transformations/geometric.py | 17 +- test/datapipes/readers/test_mesh_readers.py | 17 +- .../transforms/test_mesh_to_domain_mesh.py | 39 ++- .../transforms/test_nested_fields.py | 4 +- test/examples/test_globe_measures.py | 62 ++++ test/mesh/calculus/test_measure.py | 60 ++-- test/mesh/calculus/test_point_measures.py | 324 ++++++++++++++++++ .../io/io_zarr/test_reader_integration.py | 31 ++ test/models/globe/test_measure.py | 10 +- 25 files changed, 1040 insertions(+), 280 deletions(-) create mode 100644 test/examples/test_globe_measures.py create mode 100644 test/mesh/calculus/test_point_measures.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a51fd82755..c6ec81fd71 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,8 +26,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 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 - effective measure (area times any composed measure weights) on the interior - under the reserved `point_data` key `TARGET_QUADRATURE_MEASURE_KEY`, so + 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 @@ -36,6 +36,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### 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. + - `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 diff --git a/docs/api/mesh/calculus.rst b/docs/api/mesh/calculus.rst index f8497c4f75..0a48170c43 100644 --- a/docs/api/mesh/calculus.rst +++ b/docs/api/mesh/calculus.rst @@ -50,6 +50,124 @@ 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 the complete measure associated +with each cell or point sample. It has shape ``(n_cells,)`` in ``cell_data`` or +``(n_points,)`` in ``point_data``. It is never a dimensionless correction that +must still be multiplied by a geometric area or volume. + +``cell_measures(mesh)`` returns the explicit cell measures, or the geometric +simplex measures when the field is absent. ``point_measures(mesh)`` requires +explicit point measures: it never changes its interpretation based on whether +``mesh.cells`` is empty. An ordinary sum is counting measure; it can also be +represented explicitly by installing ones with ``dimension=0``. + +.. code:: python + + from physicsnemo.mesh.calculus import ( + cell_measures, + point_measures, + scale_measures, + set_point_measures, + ) + + # A sampling stage retains k of N cells. Prior corrections are preserved. + sampled = mesh.slice_cells(indices) + scale_measures(sampled, mesh.n_cells / sampled.n_cells) + + # Transfer complete measures when cells become independent point samples. + queries = Mesh(points=sampled.cell_centroids) + set_point_measures( + queries, cell_measures(sampled), dimension=sampled.n_manifold_dims + ) + integral = queries.integrate_samples(predictions) + +There are two distinct integration operations: + +* ``mesh.integrate(field, data_source="cells")`` integrates piecewise-constant + cell values over the cells. ``data_source="points"`` integrates a + piecewise-linear vertex field over those same cells. Both use effective + **cell** measures. Vertex fields keep the existing rule that a NaN vertex + invalidates its incident cell's contribution when ``nan_policy="omit"``. +* ``mesh.integrate_samples(field)`` sums independent point samples times explicit + **point** measures, with no dependence on connectivity. NaN omission applies + independently to each sample contribution. Missing point measures raise an + error; there is no implicit counting or geometric fallback. + +``lumped_point_measures(mesh)`` explicitly constructs vertex quadrature by +sharing each cell's effective measure equally among its vertices. For finite +nodal fields it reproduces piecewise-linear integration. Calling it does not +modify the mesh or change either integration rule. + +Measure lifecycle +~~~~~~~~~~~~~~~~~ + +All storage and reweighting helpers live in ``physicsnemo.mesh.calculus.measure``. +``set_cell_measures`` and ``set_point_measures`` assign complete measures; +``scale_measures`` multiplies existing measures by a scalar or per-entity factor. +Sampling uses the latter with inverse inclusion probabilities. Raw slicing is a +restriction and does not apply a sampling correction. Serialization and device +transfers preserve the fields. + +``mesh.to_point_cloud(point_source="cell_centroids")`` and ``mesh.to_dual_graph`` +transfer each cell's complete measure to its centroid, including the represented +dimension. Dual-graph edges retain their own geometric length measure; the +original surface or volume measure is associated with the graph's points. + +Point measures carry a scalar ``_point_measure_dimension`` in their mesh's +``global_data``. ``set_point_measures`` writes this metadata: 0 for counting, 1 +for length, 2 for area, and 3 for volume. It describes the represented measure, +not the point cloud's topological dimension. Merging point quadrature requires +matching measure dimensions and preserves this scalar metadata. + +Rigid transformations preserve measures. Uniform scaling by ``s`` multiplies +measures of dimension ``d`` by ``abs(s)**d``. Cell geometry changes preserve the +ratio of represented to geometric measure. Subdivision transfers that ratio to +children; linear subdivision therefore conserves each parent's total measure. +A nonzero measure on a geometrically degenerate cell needs explicit replacement +measures when its geometry changes. + +For points, full-dimensional measures also support general square linear maps +through their absolute determinant. Anisotropic transformations of embedded +surface/curve quadrature require support geometry that a point cloud does not +contain, and are rejected. Transform the source cells before creating those +samples, or explicitly retain reference measures with +``mesh.with_points(new_points, preserve_measures=True)``. This policy also serves +coordinate normalization where measures deliberately remain in reference units. + +Generic field interpolation excludes effective measures. Subdivision of explicit +point quadrature and remeshing of explicit cell/point quadrature require an +explicit conservative transfer or replacement measures; ordinary field +interpolation cannot provide one. + +This pre-release API replaces ``_measure_weights`` and the datapipes-specific +target quadrature key. Files using the old cell multiplier must be regenerated +or converted to complete measures before use. Existing geometric areas remain +geometric: producers must not override area caches to store represented areas. + +For a legacy mesh, convert the stored multiplier once and remove its old key: + +.. code:: python + + from physicsnemo.mesh.calculus import set_cell_measures, set_point_measures + + if "_measure_weights" in mesh.cell_data: + weights = mesh.cell_data.pop("_measure_weights") + set_cell_measures(mesh, mesh.cell_areas * weights) + + # Legacy centroid query measures already contain the geometric contribution. + # Use the dimension of the source cells (2 here for a surface), not the + # zero-dimensional topology of the query cloud. + if "_target_quadrature_measure" in mesh.point_data: + measures = mesh.point_data.pop("_target_quadrature_measure") + set_point_measures(mesh, measures, dimension=2) + +Update field mappings to ``cell_data._effective_measure`` or +``point_data._effective_measure`` as appropriate. Do not multiply these complete +measures by geometric areas a second time. + API Reference ------------- diff --git a/examples/cfd/external_aerodynamics/globe/drivaer/dataset.py b/examples/cfd/external_aerodynamics/globe/drivaer/dataset.py index cfdd5849c3..4f236f7ca1 100644 --- a/examples/cfd/external_aerodynamics/globe/drivaer/dataset.py +++ b/examples/cfd/external_aerodynamics/globe/drivaer/dataset.py @@ -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 diff --git a/physicsnemo/datapipes/readers/mesh.py b/physicsnemo/datapipes/readers/mesh.py index cc0d5eba4a..2a7f3892c6 100644 --- a/physicsnemo/datapipes/readers/mesh.py +++ b/physicsnemo/datapipes/readers/mesh.py @@ -35,7 +35,7 @@ from physicsnemo.datapipes._rng import spawn_generator from physicsnemo.datapipes.registry import register from physicsnemo.mesh import DomainMesh, Mesh -from physicsnemo.mesh.calculus.measure import compose_measure_weights +from physicsnemo.mesh.calculus.measure import EFFECTIVE_MEASURE_KEY, scale_measures logger = logging.getLogger(__name__) @@ -60,7 +60,7 @@ def _subsample_mesh_points( falls back to ``slice_points``. Unlike :func:`_subsample_mesh_cells`, this does NOT maintain - measure weights: dropping points removes cells implicitly, with + cell-measure corrections: dropping points removes cells implicitly, with no per-cell inclusion probability to invert. Prefer cell subsampling when downstream code integrates over the mesh. """ @@ -73,14 +73,18 @@ def _subsample_mesh_points( device=mesh.points.device, ) if mesh.n_cells == 0: - return Mesh( + result = Mesh( points=mesh.points[indices], cells=mesh.cells, point_data=mesh.point_data[indices], cell_data=mesh.cell_data, global_data=mesh.global_data, ) - return mesh.slice_points(indices) + else: + result = mesh.slice_points(indices) + if EFFECTIVE_MEASURE_KEY in result.point_data: + scale_measures(result, mesh.n_points / n_points, association="points") + return result def _subsample_mesh_cells( @@ -97,8 +101,8 @@ def _subsample_mesh_cells( Preserves the mesh's integration measure: every cell's inclusion probability is exactly ``k/N``, and the retained cells' measure - weights (see :mod:`physicsnemo.mesh.calculus.measure`) are multiplied by - ``N/k``, composing with any weights from earlier sampling stages. + (see :mod:`physicsnemo.mesh.calculus.measure`) is multiplied by + ``N/k``, composing with corrections from earlier sampling stages. Consumers of the effective cell measure (see :mod:`physicsnemo.mesh.calculus.measure`) then see an unbiased estimate of the full-mesh measure rather than the ~``k/N`` retained fraction. @@ -125,7 +129,7 @@ def _subsample_mesh_cells( ### Compose the Horvitz-Thompson weight for this sampling stage. ### slice_cells/slice_points returned fresh TensorDicts, so the ### in-place update cannot leak into the memmap-backed source. - compose_measure_weights(mesh, n_total / n_cells) + scale_measures(mesh, n_total / n_cells) return mesh @@ -146,7 +150,7 @@ def _zarr_mesh_subsampled( """Partial-read a zarr mesh group: fetch only the subsample window. Reproduces :func:`_subsample_mesh` semantics (cyclic contiguous blocks, - vertex compaction, Horvitz-Thompson measure weights) while reading only + vertex compaction, Horvitz-Thompson measure corrections) while reading only the selected rows from the store instead of materializing the full mesh. """ from physicsnemo.mesh.io import io_zarr as _ioz @@ -174,7 +178,7 @@ def _zarr_mesh_subsampled( ), global_data=_ioz._read_tree(group, "global_data"), ) - compose_measure_weights(mesh, total_cells / n_cells) + scale_measures(mesh, total_cells / n_cells) if n_points is not None: mesh = _subsample_mesh_points(mesh, n_points, generator=generator) return mesh @@ -182,7 +186,7 @@ def _zarr_mesh_subsampled( if total_cells == 0 and n_points is not None and total_points > n_points: indices = _cyclic_block_indices(total_points, n_points, generator=generator) runs = _indices_to_runs(indices) - return Mesh( + mesh = Mesh( points=_ioz._read_rows(group["points"], runs), point_data=_ioz._read_tree( group, "point_data", leaf_reader=lambda a: _ioz._read_rows(a, runs) @@ -190,6 +194,9 @@ def _zarr_mesh_subsampled( cell_data=_ioz._read_tree(group, "cell_data"), global_data=_ioz._read_tree(group, "global_data"), ) + if EFFECTIVE_MEASURE_KEY in mesh.point_data: + scale_measures(mesh, total_points / n_points, association="points") + return mesh # No subsampling applies (small mesh, or unsupported combination): # eager full read keeps semantics identical to the memmap path. @@ -263,7 +270,7 @@ def __init__( choice for triangulated surface meshes where downstream transforms depend on cells (e.g. surface normals, cell centroids, cell_data fields). Records the inverse inclusion - probability as measure weights, preserving the integration + probability to effective measures, preserving the integration measure (see :mod:`physicsnemo.mesh.calculus.measure`). Applied before ``subsample_n_points`` when both are set. """ @@ -449,7 +456,7 @@ def __init__( sequential I/O, then compacts unreferenced vertices. Preserves cell topology and is the correct choice when downstream transforms depend on cells. Records the - inverse inclusion probability as measure weights, preserving + inverse inclusion probability to effective measures, preserving the integration measure (see :mod:`physicsnemo.mesh.calculus.measure`). Applied before diff --git a/physicsnemo/datapipes/transforms/mesh/__init__.py b/physicsnemo/datapipes/transforms/mesh/__init__.py index b032bb46c9..0fd65ca02a 100644 --- a/physicsnemo/datapipes/transforms/mesh/__init__.py +++ b/physicsnemo/datapipes/transforms/mesh/__init__.py @@ -29,7 +29,6 @@ ) from physicsnemo.datapipes.transforms.mesh.base import MeshTransform from physicsnemo.datapipes.transforms.mesh.transforms import ( - TARGET_QUADRATURE_MEASURE_KEY, CenterMesh, ComputeCellCentroids, ComputeSurfaceNormals, @@ -60,7 +59,6 @@ "SetGlobalField", "NormalizeMeshFields", "MeshToDomainMesh", - "TARGET_QUADRATURE_MEASURE_KEY", "MeshToTensorDict", "RestructureTensorDict", "RandomScaleMesh", diff --git a/physicsnemo/datapipes/transforms/mesh/transforms.py b/physicsnemo/datapipes/transforms/mesh/transforms.py index a729da0bec..5787cf8a96 100644 --- a/physicsnemo/datapipes/transforms/mesh/transforms.py +++ b/physicsnemo/datapipes/transforms/mesh/transforms.py @@ -45,16 +45,16 @@ Mesh, MeshFieldAssociation, ) -from physicsnemo.mesh.calculus.measure import cell_measures, compose_measure_weights +from physicsnemo.mesh.calculus.measure import ( + EFFECTIVE_MEASURE_KEY, + cell_measures, + point_measure_dimension, + point_measures, + scale_measures, + set_point_measures, +) from physicsnemo.nn.functional import weighted_multinomial -### Reserved ``point_data`` key carrying the effective quadrature measure of -### the cell-centroid query points created by :class:`MeshToDomainMesh`. -### Distinct from ``MEASURE_WEIGHTS_KEY``: that is a dimensionless factor on -### source cells, while this is the full geometric measure (area * weight), -### aligned one-for-one with the interior points. -TARGET_QUADRATURE_MEASURE_KEY: str = "_target_quadrature_measure" - @register() class ScaleMesh(MeshTransform): @@ -275,10 +275,10 @@ def _compact_points(mesh: Mesh) -> Mesh: class SubsampleMesh(MeshTransform): r"""Subsample a mesh to a fixed number of cells and/or points. - Cell subsampling preserves the integration measure by recording - each stage's inverse inclusion probability into the mesh's measure - weights (see :mod:`physicsnemo.mesh.calculus.measure`); point - subsampling does not maintain weights. + Sampling multiplies explicit effective measures by the stage's inverse + inclusion probability. Cell measures default to geometric measures; + point measures remain absent unless supplied explicitly. Point sampling + does not correct cell measures for cells removed indirectly. """ def __init__( @@ -324,16 +324,19 @@ def __call__(self, mesh: Mesh) -> Mesh: if self.compact: mesh = _compact_points(mesh) ### Compose this stage's inverse inclusion probability into the - ### mesh's measure weights. + ### mesh's effective measures. ### `_random_indices` is exact below the large-population threshold ### and uses the near-uniform Poisson-gap approximation above it. - compose_measure_weights(mesh, n_before / self.n_cells) + scale_measures(mesh, n_before / self.n_cells) if self.n_points is not None and mesh.n_points > self.n_points: indices = self._random_indices( mesh.n_points, self.n_points, mesh.points.device ) + n_before = mesh.n_points mesh = mesh.slice_points(indices) + if EFFECTIVE_MEASURE_KEY in mesh.point_data: + scale_measures(mesh, n_before / self.n_points, association="points") return mesh @@ -1075,7 +1078,7 @@ class MeshToDomainMesh(MeshTransform): If ``None`` (and ``point_data_targets`` is also ``None``), no user targets are placed on the interior. Centroid mode still records each source cell's effective measure under - :data:`TARGET_QUADRATURE_MEASURE_KEY`, so integrals over the query + :data:`~physicsnemo.mesh.calculus.measure.EFFECTIVE_MEASURE_KEY`, so integrals over the query points remain possible after the cells are gone. point_data_targets : list[str] or None, default ``None`` Names of vertex-centered fields on the input mesh to use as prediction @@ -1158,27 +1161,16 @@ def __init__( ### ``select`` / ``exclude`` below accept the parsed tuple keys. self._cell_data_targets: list[NestedKey] = as_nested_keys(cell_data_targets) self._point_data_targets: list[NestedKey] = as_nested_keys(point_data_targets) - if TARGET_QUADRATURE_MEASURE_KEY in ( - self._cell_data_targets + self._point_data_targets - ): + reserved = as_nested_key(EFFECTIVE_MEASURE_KEY) + if reserved in self._cell_data_targets or reserved in self._point_data_targets: raise ValueError( - f"{TARGET_QUADRATURE_MEASURE_KEY!r} is reserved for the query " - "measure and cannot be configured as a target." + f"{EFFECTIVE_MEASURE_KEY!r} is reserved for effective " + "measure bookkeeping and cannot be configured as a user target." ) self._interior_points = interior_points self._boundary_name = boundary_name def __call__(self, mesh: Mesh) -> DomainMesh: # type: ignore[override] - for association, data in ( - ("point_data", mesh.point_data), - ("cell_data", mesh.cell_data), - ): - if TARGET_QUADRATURE_MEASURE_KEY in data: - raise ValueError( - f"Input mesh {association} already contains reserved key " - f"{TARGET_QUADRATURE_MEASURE_KEY!r}; rename the field " - "before MeshToDomainMesh." - ) ### v1 supports two diagonal corners: ### (cell_data_targets, interior_points='cell_centroids') ### (point_data_targets, interior_points='vertices') @@ -1212,20 +1204,23 @@ def __call__(self, mesh: Mesh) -> DomainMesh: # type: ignore[override] def _call_cell_centroids(self, mesh: Mesh) -> DomainMesh: ### Build the interior as a point cloud at cell centroids, with target - ### cell_data fields moved into interior.point_data. The source cells - ### do not exist on the interior, so record their effective measure - ### (area * any composed measure weights) beside the centroids now. + ### cell_data fields moved into interior.point_data. The original + ### cells disappear at this boundary, so materialize their effective + ### measure beside the centroid queries while cell geometry and any + ### sampling corrections are still available. require_keys(mesh.cell_data, self._cell_data_targets, what="Target field") interior_point_data = ( mesh.cell_data.select(*self._cell_data_targets) if self._cell_data_targets else TensorDict({}, batch_size=[mesh.n_cells]) ) - interior_point_data[TARGET_QUADRATURE_MEASURE_KEY] = cell_measures(mesh) interior = Mesh( points=mesh.cell_centroids, point_data=interior_point_data, ) + set_point_measures( + interior, cell_measures(mesh), dimension=mesh.n_manifold_dims + ) ### Build the boundary by stripping target fields from cell_data. boundary_cell_data = ( exclude_keys(mesh.cell_data, self._cell_data_targets) @@ -1252,6 +1247,12 @@ def _call_vertices(self, mesh: Mesh) -> DomainMesh: points=mesh.points, point_data=interior_point_data, ) + if EFFECTIVE_MEASURE_KEY in mesh.point_data: + set_point_measures( + interior, + point_measures(mesh), + dimension=int(point_measure_dimension(mesh)), + ) boundary_point_data = ( exclude_keys(mesh.point_data, self._point_data_targets) if self._point_data_targets diff --git a/physicsnemo/mesh/calculus/__init__.py b/physicsnemo/mesh/calculus/__init__.py index 48e5d46a64..6f9d6817cd 100644 --- a/physicsnemo/mesh/calculus/__init__.py +++ b/physicsnemo/mesh/calculus/__init__.py @@ -59,13 +59,18 @@ integrate_flux, integrate_moment, integrate_point_data, + integrate_samples, ) from physicsnemo.mesh.calculus.laplacian import ( compute_laplacian_points_dec, ) from physicsnemo.mesh.calculus.measure import ( - MEASURE_WEIGHTS_KEY, - cell_measure_weights, + EFFECTIVE_MEASURE_KEY, + POINT_MEASURE_DIMENSION_KEY, cell_measures, - compose_measure_weights, + lumped_point_measures, + point_measures, + scale_measures, + set_cell_measures, + set_point_measures, ) diff --git a/physicsnemo/mesh/calculus/integration.py b/physicsnemo/mesh/calculus/integration.py index 3a1b75b9f3..2f96fd8b51 100644 --- a/physicsnemo/mesh/calculus/integration.py +++ b/physicsnemo/mesh/calculus/integration.py @@ -55,7 +55,7 @@ from jaxtyping import Float from physicsnemo.core.warnings import LegacyFeatureWarning -from physicsnemo.mesh.calculus.measure import cell_measures +from physicsnemo.mesh.calculus.measure import cell_measures, point_measures if TYPE_CHECKING: from physicsnemo.mesh.mesh import Mesh @@ -320,6 +320,28 @@ def integrate( raise ValueError(f"Invalid {data_source=!r}. Must be 'cells' or 'points'.") +def integrate_samples( + mesh: "Mesh", + field: str | tuple[str, ...] | torch.Tensor, + *, + nan_policy: NanPolicy = "omit", +) -> torch.Tensor: + """Integrate independent point samples using explicit point measures. + + This is a quadrature sum, regardless of whether the mesh has cells. It + does not interpolate values through cells or infer a point measure. + Missing measures raise KeyError. Use an ordinary tensor sum for counting + measure, or set_point_measures with ones and dimension=0 explicitly. + For a piecewise-linear vertex field over connected cells, use integrate + with data_source="points" instead; that retains cell-based NaN handling. + """ + values = _resolve_field(mesh, field, "points") + if values.ndim == 0 or values.shape[0] != mesh.n_points: + raise ValueError("Point sample leading dimension must equal n_points") + measures = point_measures(mesh).reshape(-1, *([1] * (values.ndim - 1))) + return _sum_with_nan_policy(values * measures, dim=0, nan_policy=nan_policy) + + def integrate_cell_data( mesh: "Mesh", field: Float[torch.Tensor, "n_cells ..."], diff --git a/physicsnemo/mesh/calculus/measure.py b/physicsnemo/mesh/calculus/measure.py index 59bd6c5f8a..14cf1830d1 100644 --- a/physicsnemo/mesh/calculus/measure.py +++ b/physicsnemo/mesh/calculus/measure.py @@ -14,150 +14,228 @@ # See the License for the specific language governing permissions and # limitations under the License. -r"""Discrete integration measure for meshes. - -A mesh discretizes a manifold, and every integral over it is a -measure-weighted sum: each cell contributes its field value times the -measure that cell represents. For a mesh that represents exactly the -geometry it stores, that measure is the geometric simplex measure, -:attr:`~physicsnemo.mesh.Mesh.cell_areas`. But the two can diverge -- a -cell may stand for more, less, or other than itself. This module defines -the resulting contract: the effective **cell measure** - -.. math:: - \mu_c = |\sigma_c| \, w_c, - -where :math:`|\sigma_c|` is the geometric measure of cell :math:`c` and -:math:`w_c` is a dimensionless per-cell **measure weight** -- the ratio of -represented to geometric measure -- defaulting to one. - -Non-unit weights arise whenever representation and geometry part ways: -cells standing in for symmetric or periodic images of themselves, -fractional weights that de-duplicate overlapping patches, corrections for -curved geometry that a straight simplex under-resolves, or coarse cells -representing a partition of finer ones. The canonical source in this -package is random cell subsampling: a cell retained with inclusion -probability :math:`\pi_c` statistically represents :math:`1/\pi_c` cells -of the original mesh, so the subsampling stage records the -Horvitz-Thompson weight :math:`w_c = 1/\pi_c` (exactly ``N/k`` for a -uniform ``k``-of-``N`` subsample). Integrals computed with the effective -measure are then unbiased estimates of the full-mesh integrals; computed -with the bare geometric measure they shrink by ``~k/N``, compounded once -per stage in any computation that chains measure-weighted sums. - -The contract has three parts: - -- Measure weights are stored in ``cell_data`` under the reserved key - :data:`MEASURE_WEIGHTS_KEY`. Living in ``cell_data``, they survive - cell slicing, serialization, device transfer, and rigid transforms - automatically; being dimensionless, they are also invariant under - geometric rescaling (``cell_areas`` alone picks up the appropriate power - of length). The underscore prefix marks the field as bookkeeping: - feature-selection code that consumes ``cell_data`` wholesale should - exclude it. -- Producers record their weight contribution via - :func:`compose_measure_weights`; successive contributions compose - multiplicatively, so e.g. chained subsampling stays exact. Point - subsampling on meshes with cells does **not** maintain weights (cells - dropped implicitly have no per-cell inclusion probability). -- Integral consumers (:meth:`Mesh.integrate`, :meth:`Mesh.integrate_flux`, - :func:`~physicsnemo.mesh.calculus.integration.integrate_moment`) read - :func:`cell_measures`. Meshes without recorded weights pass through - with the bare geometric measure, bit-identically. +r"""Complete integration measures for cells and point samples. + +``_effective_measure`` always stores the complete contribution to an integral, +never a multiplier that must still be multiplied by geometry. Cells without an +explicit measure use their geometric simplex measure. Point measures are always +explicit, independent of connectivity; counting measure is an explicit choice. + +Sampling producers multiply effective measures by their inverse inclusion +probability. Cell-to-centroid conversion transfers those measures to points. +All such producers use the helpers here rather than inventing storage keys. + +Point quadrature also records its represented dimension in ``global_data``: +zero for counting, one for length, two for area, etc. This describes the measure, +not the topology of the point cloud. It suffices for uniform scaling; arbitrary +changes of point positions require support geometry or an explicit decision to +preserve the measure. Cell measures follow the ratio of new to old geometric +measure when geometry changes, preserving their represented/geometric ratio. """ -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Literal import torch -from jaxtyping import Float if TYPE_CHECKING: from physicsnemo.mesh.mesh import Mesh -### Reserved `cell_data` key holding dimensionless per-cell measure weights -### (ratios of represented to geometric measure). See the module docstring -### for the contract. -MEASURE_WEIGHTS_KEY: str = "_measure_weights" +EFFECTIVE_MEASURE_KEY = "_effective_measure" +POINT_MEASURE_DIMENSION_KEY = "_point_measure_dimension" -def _validate_measure_weight_shape( - mesh: "Mesh", weights: torch.Tensor -) -> Float[torch.Tensor, " n_cells"]: - """Enforce the reserved field's one-scalar-per-cell contract. +def _validate_measures( + mesh: "Mesh", values: torch.Tensor, association: str +) -> torch.Tensor: + n = mesh.n_cells if association == "cells" else mesh.n_points + if not isinstance(values, torch.Tensor) or values.shape != (n,): + shape = ( + tuple(values.shape) + if isinstance(values, torch.Tensor) + else type(values).__name__ + ) + raise ValueError( + f"{EFFECTIVE_MEASURE_KEY!r} must have shape {(n,)} (one scalar per {association}), got {shape}" + ) + if not values.is_floating_point(): + raise TypeError("Effective measures must be floating-point tensors") + if values.device != mesh.points.device: + raise ValueError( + "Effective measures and mesh points must be on the same device" + ) + return values + - ``cell_data``'s batch dimension only checks the leading axis, so a - ``(n_cells, 1)`` tensor is storable but would broadcast wrongly against - ``cell_areas``. +def cell_measures(mesh: "Mesh") -> torch.Tensor: + """Return complete per-cell measures, falling back to geometric measures. + + The result has shape ``(n_cells,)``. The fallback does not materialize a + stored field, so unweighted meshes retain lazy geometric computation. """ - expected = (mesh.n_cells,) - if weights.shape != expected: + if "_measure_weights" in mesh.cell_data: raise ValueError( - f"{MEASURE_WEIGHTS_KEY!r} must have shape {expected} " - f"(one scalar per cell), got {tuple(weights.shape)}" + "Legacy _measure_weights must be converted to _effective_measure = " + "cell_areas * weights before use; the new field stores complete measures" ) - return weights + values = mesh.cell_data.get(EFFECTIVE_MEASURE_KEY, None) + return ( + mesh.cell_areas if values is None else _validate_measures(mesh, values, "cells") + ) -def cell_measure_weights(mesh: "Mesh") -> Float[torch.Tensor, " n_cells"]: - r"""Per-cell measure weights of *mesh* (ones when none are recorded). +def point_measures(mesh: "Mesh") -> torch.Tensor: + """Return explicit per-point measures of shape ``(n_points,)``. - Returns - ------- - torch.Tensor - Dimensionless weights of shape ``(n_cells,)``. If no weights have - been recorded, returns ones: every cell represents exactly itself. + Raises ``KeyError`` when absent, even if cells exist. Use an ordinary sum + for counting measure, or explicitly install ones with dimension zero. + Use :func:`lumped_point_measures` to construct nodal quadrature from cells. """ - weights = mesh.cell_data.get(MEASURE_WEIGHTS_KEY, None) - if weights is None: - return torch.ones( - mesh.n_cells, dtype=mesh.points.dtype, device=mesh.points.device + values = mesh.point_data.get(EFFECTIVE_MEASURE_KEY, None) + if values is None: + raise KeyError( + "Point quadrature requires explicit effective measures; use set_point_measures or an ordinary sum for counting measure" ) - return _validate_measure_weight_shape(mesh, weights) + return _validate_measures(mesh, values, "points") + +def set_cell_measures(mesh: "Mesh", values: torch.Tensor) -> None: + """Set complete per-cell measures in place, without changing geometry.""" + mesh.cell_data[EFFECTIVE_MEASURE_KEY] = _validate_measures(mesh, values, "cells") -def cell_measures(mesh: "Mesh") -> Float[torch.Tensor, " n_cells"]: - r"""Effective per-cell integration measure: ``cell_areas * measure_weights``. - This is what integral consumers should weight by. Skips the - multiplication when no measure weights are recorded, so meshes without - weights pay nothing and results are bit-identical to the bare geometric - measure. +def set_point_measures(mesh: "Mesh", values: torch.Tensor, *, dimension: int) -> None: + """Set complete point measures and their represented dimension in place. - Returns - ------- - torch.Tensor - Effective measure of shape ``(n_cells,)``. + ``dimension`` is the power of length: 0 for counting, 1 for length, 2 for + area, 3 for volume. It never depends on whether this mesh has cells. """ - cell_areas = mesh.cell_areas - weights = mesh.cell_data.get(MEASURE_WEIGHTS_KEY, None) - if weights is None: - return cell_areas - return cell_areas * _validate_measure_weight_shape(mesh, weights) - - -def compose_measure_weights(mesh: "Mesh", factor: float | torch.Tensor) -> None: - r"""Multiply *mesh*'s measure weights by *factor*, in place. - - Called by each producer with its contribution to the represented-to- - geometric measure ratio -- a sampling stage, for example, passes its - inverse inclusion probability (``n_cells_before / n_cells_after`` for - a uniform sample). Contributions compose multiplicatively: a stage - keeping ``k1`` of ``N`` cells followed by one keeping ``k2`` of ``k1`` - yields exactly ``N/k2``. - - Parameters - ---------- - mesh : Mesh - Mesh to update. Its ``cell_data`` is modified in place. - factor : float or torch.Tensor - This producer's weight contribution; a scalar or a per-cell tensor - broadcast against the existing weights. + if ( + isinstance(dimension, bool) + or not isinstance(dimension, int) + or not 0 <= dimension <= mesh.n_spatial_dims + ): + raise ValueError( + f"Measure dimension must be an integer in [0, {mesh.n_spatial_dims}]" + ) + values = _validate_measures(mesh, values, "points") + mesh.point_data[EFFECTIVE_MEASURE_KEY] = values + mesh.global_data[POINT_MEASURE_DIMENSION_KEY] = torch.tensor( + dimension, device=mesh.points.device + ) + + +def point_measure_dimension(mesh: "Mesh") -> torch.Tensor: + """Read the scalar dimension metadata required to transform point measures.""" + dimension = mesh.global_data.get(POINT_MEASURE_DIMENSION_KEY, None) + if ( + dimension is None + or not isinstance(dimension, torch.Tensor) + or dimension.ndim != 0 + ): + raise ValueError( + "Point measure dimension is missing or not scalar; use set_point_measures" + ) + if not torch.compiler.is_compiling() and not bool( + (dimension >= 0) + & (dimension <= mesh.n_spatial_dims) + & (dimension == dimension.round()) + ): + raise ValueError( + "Point measure dimension must be a nonnegative integer no larger than the spatial dimension" + ) + return dimension + + +def scale_measures( + mesh: "Mesh", + factor: float | torch.Tensor, + *, + association: Literal["cells", "points"] = "cells", +) -> None: + """Multiply effective measures in place by a scalar or per-entity factor. + + For example, a sampling stage keeping k of N cells uses N/k; subsequent + stages multiply the measures already recorded. Point measures must exist + before reweighting: this function never invents a point base measure. + """ + if association not in ("cells", "points"): + raise ValueError("association must be 'cells' or 'points'") + values = cell_measures(mesh) if association == "cells" else point_measures(mesh) + if ( + isinstance(factor, torch.Tensor) + and factor.ndim != 0 + and factor.shape != values.shape + ): + raise ValueError( + f"Measure factor must be a scalar or have shape {tuple(values.shape)}, got {tuple(factor.shape)}" + ) + getattr(mesh, "cell_data" if association == "cells" else "point_data")[ + EFFECTIVE_MEASURE_KEY + ] = values * factor + + +def lumped_point_measures(mesh: "Mesh") -> torch.Tensor: + """Construct P1 nodal quadrature by dividing each cell's measure equally. + + Each vertex receives the sum of its incident cells' contributions. This is + an explicit construction, not the default point measure. For finite nodal + values it reproduces the piecewise-linear integral. Cell-based NaN omission + still requires :func:`integrate`, rather than a sum of nodal contributions. """ - if isinstance(factor, torch.Tensor) and factor.ndim != 0: - expected = (mesh.n_cells,) - if factor.shape != expected: + measures = cell_measures(mesh) + result = measures.new_zeros(mesh.n_points) + shares = (measures / mesh.cells.shape[1]).unsqueeze(-1).expand_as(mesh.cells) + return result.scatter_add(0, mesh.cells.flatten(), shares.flatten()) + + +def _transfer_cell_measures( + source: "Mesh", result: "Mesh", parents: torch.Tensor | None = None +) -> None: + """Transfer explicit cell measures through geometric changes/subdivision.""" + if EFFECTIVE_MEASURE_KEY not in source.cell_data: + return + measures = cell_measures(source) + old_areas = source.cell_areas + if parents is not None: + measures, old_areas = measures[parents], old_areas[parents] + if not torch.compiler.is_compiling() and bool( + ((old_areas == 0) & (measures != 0)).any() + ): + raise ValueError( + "Cannot transform nonzero effective measure on a zero-measure cell; supply replacement measures" + ) + denominator = torch.where(old_areas != 0, old_areas, torch.ones_like(old_areas)) + set_cell_measures(result, measures * (result.cell_areas / denominator)) + + +def _require_preserved_point_measures(mesh: "Mesh") -> None: + """Reject unknown geometric changes of dimensional point quadrature.""" + if EFFECTIVE_MEASURE_KEY in mesh.point_data and bool( + point_measure_dimension(mesh) != 0 + ): + raise ValueError( + "Changing quadrature point geometry requires replacement measures or an explicit preservation policy; use with_points(..., preserve_measures=True) to retain reference measures" + ) + + +def _transform_point_measures( + source: "Mesh", result: "Mesh", matrix: torch.Tensor +) -> None: + """Transform counting, full-dimensional, or similarity-mapped point measures.""" + if EFFECTIVE_MEASURE_KEY not in source.point_data: + return + dimension = point_measure_dimension(source) + if bool(dimension == 0): + return + if matrix.shape[0] == matrix.shape[1] and bool(dimension == source.n_spatial_dims): + factor = matrix.det().abs() + else: + from physicsnemo.mesh.transformations.geometric import _is_similarity_transform + + if not _is_similarity_transform(matrix): raise ValueError( - f"measure-weight factor must be a scalar or have shape " - f"{expected}, got {tuple(factor.shape)}" + "Anisotropic transformation of point measures requires support geometry; transform the source cells before constructing quadrature, or explicitly preserve reference measures" ) - mesh.cell_data[MEASURE_WEIGHTS_KEY] = cell_measure_weights(mesh) * factor + length_scale = (matrix.T @ matrix).diagonal().mean().clamp_min(0).sqrt() + factor = length_scale**dimension + result.point_data[EFFECTIVE_MEASURE_KEY] = point_measures(source) * factor diff --git a/physicsnemo/mesh/mesh.py b/physicsnemo/mesh/mesh.py index 6fa6f6391e..9a4c6c4c29 100644 --- a/physicsnemo/mesh/mesh.py +++ b/physicsnemo/mesh/mesh.py @@ -45,6 +45,7 @@ integrate, integrate_flux, integrate_moment, + integrate_samples, ) from physicsnemo.mesh.geometry._cell_areas import compute_cell_areas from physicsnemo.mesh.geometry._cell_normals import compute_cell_normals @@ -1401,8 +1402,29 @@ def merge( cell_index_offsets = cumsum_n_points.roll(1) cell_index_offsets[0] = 0 + from physicsnemo.mesh.calculus.measure import ( + EFFECTIVE_MEASURE_KEY, + POINT_MEASURE_DIMENSION_KEY, + point_measure_dimension, + ) + + point_dimension = None + if EFFECTIVE_MEASURE_KEY in meshes[0].point_data: + point_dimension = point_measure_dimension(meshes[0]) + if any( + not torch.equal(point_measure_dimension(m), point_dimension) + for m in meshes[1:] + ): + raise ValueError( + "Cannot merge point quadrature with different measure dimensions" + ) if global_data_strategy == "stack": - global_data = TensorDict.stack([m.global_data for m in meshes]) + global_data = TensorDict.stack( + [m.global_data.exclude(POINT_MEASURE_DIMENSION_KEY) for m in meshes] + ) + if point_dimension is not None: + global_data.batch_size = [] + global_data[POINT_MEASURE_DIMENSION_KEY] = point_dimension else: raise ValueError(f"Invalid {global_data_strategy=}") @@ -1837,6 +1859,7 @@ def with_points( points: torch.Tensor, *, keep: str | tuple[str, ...] | Sequence[str | tuple[str, ...]] = "topology", + preserve_measures: builtins.bool = False, ) -> "Mesh": r"""Return a mesh with replacement point coordinates. @@ -1855,6 +1878,11 @@ def with_points( Cache keys to retain. Uses the same key semantics as :meth:`strip_caches`; defaults to the complete ``"topology"`` cache. + preserve_measures : bool, default False + Explicitly retain reference measures when replacing coordinates. + Otherwise cell measures follow geometric measure changes; dimensional + point measures require a known transformation or replacement measures. + Returns ------- Mesh @@ -1893,11 +1921,21 @@ def with_points( lambda: "with_points must preserve point indexing.", ) - return self._new_with_structure( + from physicsnemo.mesh.calculus.measure import ( + _require_preserved_point_measures, + _transfer_cell_measures, + ) + + if not preserve_measures: + _require_preserved_point_measures(self) + result = self._new_with_structure( points=points, cells=self.cells, keep=keep, ) + if not preserve_measures: + _transfer_cell_measures(self, result) + return result def with_cells( self, @@ -1962,11 +2000,15 @@ def with_cells( lambda: "with_cells must preserve simplex type.", ) - return self._new_with_structure( + result = self._new_with_structure( points=self.points, cells=cells, keep=keep, ) + from physicsnemo.mesh.calculus.measure import _transfer_cell_measures + + _transfer_cell_measures(self, result) + return result def with_data( self, @@ -2071,9 +2113,13 @@ def cell_data_to_point_data(self, overwrite_keys: builtins.bool = False) -> "Mes >>> mesh_with_point_data = mesh.cell_data_to_point_data() # doctest: +SKIP >>> # Now mesh has both cell_data["pressure"] and point_data["pressure"] """ + from physicsnemo.mesh.calculus.measure import EFFECTIVE_MEASURE_KEY + + fields = self.cell_data.exclude(EFFECTIVE_MEASURE_KEY) + # Effective measures are not interpolated; use lumped_point_measures. ### Check for key conflicts if not overwrite_keys: - src_keys = set(self.cell_data.keys(include_nested=True, leaves_only=True)) + src_keys = set(fields.keys(include_nested=True, leaves_only=True)) dst_keys = set(self.point_data.keys(include_nested=True, leaves_only=True)) conflicts = src_keys & dst_keys if conflicts: @@ -2099,7 +2145,7 @@ def cell_data_to_point_data(self, overwrite_keys: builtins.bool = False) -> "Mes self.n_cells, device=self.points.device ).repeat_interleave(n_vertices_per_cell) - converted = self.cell_data.apply( + converted = fields.apply( lambda cell_values: scatter_aggregate( src_data=cell_values[cell_indices], src_to_dst_mapping=point_indices, @@ -2150,9 +2196,12 @@ def point_data_to_cell_data(self, overwrite_keys: builtins.bool = False) -> "Mes >>> mesh_with_cell_data = mesh.point_data_to_cell_data() # doctest: +SKIP >>> # Now mesh has both point_data["temperature"] and cell_data["temperature"] """ + from physicsnemo.mesh.calculus.measure import EFFECTIVE_MEASURE_KEY + + fields = self.point_data.exclude(EFFECTIVE_MEASURE_KEY) ### Check for key conflicts if not overwrite_keys: - src_keys = set(self.point_data.keys(include_nested=True, leaves_only=True)) + src_keys = set(fields.keys(include_nested=True, leaves_only=True)) dst_keys = set(self.cell_data.keys(include_nested=True, leaves_only=True)) conflicts = src_keys & dst_keys if conflicts: @@ -2176,7 +2225,7 @@ def _mean_over_cell_vertices(point_values: torch.Tensor) -> torch.Tensor: cell_values = cell_values.to(torch.float64) return cell_values.mean(dim=1) - converted = self.point_data.apply( + converted = fields.apply( _mean_over_cell_vertices, batch_size=torch.Size([self.n_cells]), ) @@ -2412,11 +2461,12 @@ def to_dual_graph(self) -> "Mesh[1, ...]": mask = sources < targets edges = torch.stack([sources[mask], targets[mask]], dim=1) + centroids = self.to_point_cloud(point_source="cell_centroids") return Mesh( - points=self.cell_centroids, + points=centroids.points, cells=edges, - point_data=self.cell_data, - global_data=self.global_data, + point_data=centroids.point_data, + global_data=centroids.global_data, ) def to_point_cloud( @@ -2432,7 +2482,8 @@ def to_point_cloud( - ``"vertices"`` (default): Uses mesh vertices as points, preserving ``point_data``. - ``"cell_centroids"``: Uses cell centroids as points, - mapping ``cell_data`` to ``point_data``. + mapping ``cell_data`` to ``point_data``. Complete cell measures + become point measures with the source manifold's dimension. Returns ------- @@ -2457,11 +2508,20 @@ def to_point_cloud( global_data=self.global_data, ) elif point_source == "cell_centroids": - return Mesh( + from physicsnemo.mesh.calculus.measure import ( + cell_measures, + set_point_measures, + ) + + result = Mesh( points=self.cell_centroids, - point_data=self.cell_data, - global_data=self.global_data, + point_data=self.cell_data.copy(), + global_data=self.global_data.copy(), ) + set_point_measures( + result, cell_measures(self), dimension=self.n_manifold_dims + ) + return result else: raise ValueError( f"Invalid {point_source=!r}. Must be 'vertices' or 'cell_centroids'." @@ -2855,6 +2915,8 @@ def next_power_size(current_size: int, base: float) -> int: compute_point_derivatives = compute_point_derivatives + integrate_samples = integrate_samples + integrate = integrate integrate_flux = integrate_flux diff --git a/physicsnemo/mesh/remeshing/_partition.py b/physicsnemo/mesh/remeshing/_partition.py index d62dc6ac7a..aaae701854 100644 --- a/physicsnemo/mesh/remeshing/_partition.py +++ b/physicsnemo/mesh/remeshing/_partition.py @@ -132,7 +132,7 @@ def partition_cells( smooth surfaces where inter-seed spacing is small relative to the radius of curvature, this is an excellent approximation. - Every original cell is assigned to exactly one cluster, so - ``cluster_areas.sum() == mesh.cell_areas.sum()`` by construction. + ``cluster_areas.sum() == cell_measures(mesh).sum()`` by construction. - If a cluster receives no cells (possible when seeds outnumber cells or cluster heavily), its area is 0, its normal is the zero vector, and its centroid falls back to the seed position. @@ -157,7 +157,9 @@ def partition_cells( ### Read source geometry (cached on Mesh) n_dims = mesh.n_spatial_dims cell_centroids = mesh.cell_centroids # (M, D) - cell_areas = mesh.cell_areas # (M,) + from physicsnemo.mesh.calculus.measure import cell_measures + + cell_areas = cell_measures(mesh) # complete represented measures (M,) has_normals = mesh.codimension == 1 ### Assign each cell to its nearest seed via kNN search (k=1). diff --git a/physicsnemo/mesh/remeshing/_remeshing.py b/physicsnemo/mesh/remeshing/_remeshing.py index cda9270f8c..ccd8e859d3 100644 --- a/physicsnemo/mesh/remeshing/_remeshing.py +++ b/physicsnemo/mesh/remeshing/_remeshing.py @@ -299,6 +299,16 @@ def remesh( :func:`physicsnemo.nn.functional.remeshing`. These advanced parameters may change as the implementation evolves. """ + from physicsnemo.mesh.calculus.measure import EFFECTIVE_MEASURE_KEY + + if ( + EFFECTIVE_MEASURE_KEY in mesh.cell_data + or EFFECTIVE_MEASURE_KEY in mesh.point_data + ): + raise ValueError( + "Remeshing explicit quadrature requires a conservative measure transfer; supply replacement measures instead of interpolating or dropping them" + ) + if mesh.n_manifold_dims != 2 or mesh.n_spatial_dims != 3: raise NotImplementedError( "remesh only supports 2D triangle surfaces embedded in 3D. Got " diff --git a/physicsnemo/mesh/subdivision/_data.py b/physicsnemo/mesh/subdivision/_data.py index 0c7c30c0bc..d25bac9a92 100644 --- a/physicsnemo/mesh/subdivision/_data.py +++ b/physicsnemo/mesh/subdivision/_data.py @@ -66,6 +66,13 @@ def interpolate_point_data_to_edges( >>> new_data = interpolate_point_data_to_edges(point_data, edges, 3) >>> # new_data["temperature"] = [100, 200, 300, 150, 250] """ + from physicsnemo.mesh.calculus.measure import EFFECTIVE_MEASURE_KEY + + if EFFECTIVE_MEASURE_KEY in point_data: + raise ValueError( + "Subdivision of explicit point quadrature requires replacement measures; point measures cannot be interpolated as fields" + ) + if len(point_data.keys()) == 0: # No data to interpolate return TensorDict( diff --git a/physicsnemo/mesh/subdivision/butterfly.py b/physicsnemo/mesh/subdivision/butterfly.py index 6f08cba66a..bf4ff64ecd 100644 --- a/physicsnemo/mesh/subdivision/butterfly.py +++ b/physicsnemo/mesh/subdivision/butterfly.py @@ -364,10 +364,15 @@ def subdivide_butterfly(mesh: "Mesh") -> "Mesh": ) ### Create and return subdivided mesh - return Mesh( + result = Mesh( points=new_points, cells=child_cells, point_data=new_point_data, cell_data=new_cell_data, global_data=mesh.global_data, ) + + from physicsnemo.mesh.calculus.measure import _transfer_cell_measures + + _transfer_cell_measures(mesh, result, parent_indices) + return result diff --git a/physicsnemo/mesh/subdivision/linear.py b/physicsnemo/mesh/subdivision/linear.py index 480206368b..80a6172b90 100644 --- a/physicsnemo/mesh/subdivision/linear.py +++ b/physicsnemo/mesh/subdivision/linear.py @@ -127,10 +127,15 @@ def subdivide_linear(mesh: "Mesh") -> "Mesh": ) ### Create and return subdivided mesh - return Mesh( + result = Mesh( points=new_points, cells=child_cells, point_data=new_point_data, cell_data=new_cell_data, global_data=mesh.global_data, # Preserved unchanged ) + + from physicsnemo.mesh.calculus.measure import _transfer_cell_measures + + _transfer_cell_measures(mesh, result, parent_indices) + return result diff --git a/physicsnemo/mesh/subdivision/loop.py b/physicsnemo/mesh/subdivision/loop.py index 901afb5e90..4fa5ae471e 100644 --- a/physicsnemo/mesh/subdivision/loop.py +++ b/physicsnemo/mesh/subdivision/loop.py @@ -428,10 +428,15 @@ def subdivide_loop(mesh: "Mesh") -> "Mesh": ) ### Create and return subdivided mesh - return Mesh( + result = Mesh( points=new_points, cells=child_cells, point_data=new_point_data, cell_data=new_cell_data, global_data=mesh.global_data, ) + + from physicsnemo.mesh.calculus.measure import _transfer_cell_measures + + _transfer_cell_measures(mesh, result, parent_indices) + return result diff --git a/physicsnemo/mesh/transformations/geometric.py b/physicsnemo/mesh/transformations/geometric.py index 637c7ef533..d880628b23 100644 --- a/physicsnemo/mesh/transformations/geometric.py +++ b/physicsnemo/mesh/transformations/geometric.py @@ -485,6 +485,13 @@ def transform( Notes ----- + Explicit effective measures follow geometry independently of the ordinary + field-transformation flags. Cell measures follow geometric measure ratios. + Point measures use their represented dimension for similarities and the + determinant for full-dimensional square maps. Other point-measure maps + require support geometry and raise ValueError. To retain reference measures + explicitly, use ``mesh.with_points(new_points, preserve_measures=True)``. + Cache Handling: - areas: For square invertible matrices: @@ -508,7 +515,7 @@ def transform( ### Start from the cache policy for coordinate replacement: retain topology, # invalidate geometry, then opt individual transformed values back in below. - transformed_mesh = mesh.with_points(new_points) + transformed_mesh = mesh.with_points(new_points, preserve_measures=True) new_cache = transformed_mesh._cache ### Opt-in: areas and normals (only for square invertible matrices) @@ -602,6 +609,13 @@ def transform( ), ) + from physicsnemo.mesh.calculus.measure import ( + _transfer_cell_measures, + _transform_point_measures, + ) + + _transfer_cell_measures(mesh, transformed_mesh) + _transform_point_measures(mesh, transformed_mesh, matrix) return transformed_mesh @@ -648,6 +662,7 @@ def translate( new_points = mesh.points + offset translated_mesh = mesh.with_points( new_points, + preserve_measures=True, keep=( "topology", ("cell", "areas"), diff --git a/test/datapipes/readers/test_mesh_readers.py b/test/datapipes/readers/test_mesh_readers.py index 55d8b85c86..c9bbb5bebe 100644 --- a/test/datapipes/readers/test_mesh_readers.py +++ b/test/datapipes/readers/test_mesh_readers.py @@ -33,8 +33,7 @@ ) from physicsnemo.mesh import DomainMesh, Mesh from physicsnemo.mesh.calculus.measure import ( - MEASURE_WEIGHTS_KEY, - cell_measure_weights, + EFFECTIVE_MEASURE_KEY, cell_measures, ) from physicsnemo.mesh.primitives.basic import ( @@ -471,7 +470,7 @@ def test_reader_records_weight(self, tmp_path): reader.set_generator(torch.Generator().manual_seed(0)) mesh, _ = reader[0] assert mesh.n_cells == k - torch.testing.assert_close(cell_measure_weights(mesh), torch.full((k,), n / k)) + torch.testing.assert_close(cell_measures(mesh), mesh.cell_areas * (n / k)) def test_reader_noop_below_threshold(self, tmp_path): n = 8 @@ -479,7 +478,7 @@ def test_reader_noop_below_threshold(self, tmp_path): reader = MeshReader(tmp_path, pattern="*.pmsh", subsample_n_cells=20) mesh, _ = reader[0] assert mesh.n_cells == n - assert MEASURE_WEIGHTS_KEY not in mesh.cell_data.keys() + assert EFFECTIVE_MEASURE_KEY not in mesh.cell_data.keys() def test_equal_area_mesh_recovers_total_exactly(self, tmp_path): # With identical triangles, ANY cyclic block reproduces the full @@ -504,9 +503,7 @@ def test_composes_with_subsample_mesh_transform(self, tmp_path): mesh, _ = reader[0] mesh = SubsampleMesh(n_cells=k2)(mesh) assert mesh.n_cells == k2 - torch.testing.assert_close( - cell_measure_weights(mesh), torch.full((k2,), n / k2) - ) + torch.testing.assert_close(cell_measures(mesh), mesh.cell_areas * (n / k2)) def test_domain_mesh_reader_records_weights_on_boundaries(self, tmp_path): interior = Mesh(points=torch.randn(10, 3)) @@ -517,11 +514,11 @@ def test_domain_mesh_reader_records_weights_on_boundaries(self, tmp_path): loaded, _ = reader[0] assert loaded.boundaries["wall"].n_cells == 6 torch.testing.assert_close( - cell_measure_weights(loaded.boundaries["wall"]), - torch.full((6,), 24 / 6), + cell_measures(loaded.boundaries["wall"]), + loaded.boundaries["wall"].cell_areas * (24 / 6), ) ### Interior is a point cloud: no cells, no weights. - assert MEASURE_WEIGHTS_KEY not in loaded.interior.cell_data.keys() + assert EFFECTIVE_MEASURE_KEY not in loaded.interior.cell_data.keys() def test_seeded_reproducibility(self, tmp_path): ### Reader RNG is derived per-sample from (base_seed, epoch, index): diff --git a/test/datapipes/transforms/test_mesh_to_domain_mesh.py b/test/datapipes/transforms/test_mesh_to_domain_mesh.py index 147f85cd20..4c0ccde6c1 100644 --- a/test/datapipes/transforms/test_mesh_to_domain_mesh.py +++ b/test/datapipes/transforms/test_mesh_to_domain_mesh.py @@ -20,11 +20,10 @@ import torch from physicsnemo.datapipes.transforms.mesh import ( - TARGET_QUADRATURE_MEASURE_KEY, MeshToDomainMesh, ) from physicsnemo.mesh import DomainMesh, Mesh -from physicsnemo.mesh.calculus import MEASURE_WEIGHTS_KEY, cell_measures +from physicsnemo.mesh.calculus import EFFECTIVE_MEASURE_KEY, cell_measures # --------------------------------------------------------------------------- # Helpers @@ -96,7 +95,7 @@ def _unequal_triangle_mesh_3d() -> Mesh: cells=torch.tensor([[0, 1, 2], [3, 4, 5]]), cell_data={ "C_p": torch.tensor([10.0, 20.0]), - MEASURE_WEIGHTS_KEY: torch.tensor([2.0, 0.25]), + EFFECTIVE_MEASURE_KEY: torch.tensor([2.0, 0.25]), }, ) @@ -185,7 +184,11 @@ def test_targets_moved_to_interior_point_data(self): transform = MeshToDomainMesh(cell_data_targets=["C_p", "C_f"]) domain = transform(mesh) interior_keys = set(domain.interior.point_data.keys()) - assert interior_keys == {"C_p", "C_f", TARGET_QUADRATURE_MEASURE_KEY} + assert interior_keys == { + "C_p", + "C_f", + EFFECTIVE_MEASURE_KEY, + } assert torch.allclose(domain.interior.point_data["C_p"], mesh.cell_data["C_p"]) assert torch.allclose(domain.interior.point_data["C_f"], mesh.cell_data["C_f"]) @@ -195,7 +198,7 @@ def test_target_measure_is_effective_cell_measure_in_cell_order(self): boundary = domain.boundaries["vehicle"] assert torch.equal( - domain.interior.point_data[TARGET_QUADRATURE_MEASURE_KEY], + domain.interior.point_data[EFFECTIVE_MEASURE_KEY], cell_measures(mesh), ) assert torch.equal( @@ -204,12 +207,12 @@ def test_target_measure_is_effective_cell_measure_in_cell_order(self): ### Query i, target i, and measure i all refer to boundary cell i. assert torch.equal(domain.interior.points, boundary.cell_centroids) assert torch.equal( - domain.interior.point_data[TARGET_QUADRATURE_MEASURE_KEY], + domain.interior.point_data[EFFECTIVE_MEASURE_KEY], cell_measures(boundary), ) - ### The measure lives only on the interior, never on the boundary. - assert TARGET_QUADRATURE_MEASURE_KEY not in boundary.point_data - assert TARGET_QUADRATURE_MEASURE_KEY not in boundary.cell_data + ### Both associations use the same key for their own complete measure. + assert EFFECTIVE_MEASURE_KEY not in boundary.point_data + assert EFFECTIVE_MEASURE_KEY in boundary.cell_data def test_non_target_cell_data_stays_on_boundary(self): mesh = _two_triangle_mesh_3d() @@ -248,23 +251,29 @@ def test_no_targets_yields_only_target_measure(self): mesh = _two_triangle_mesh_3d() transform = MeshToDomainMesh(cell_data_targets=None) domain = transform(mesh) - assert set(domain.interior.point_data.keys()) == {TARGET_QUADRATURE_MEASURE_KEY} + assert set(domain.interior.point_data.keys()) == {EFFECTIVE_MEASURE_KEY} ### All original cell_data should still be on the boundary. boundary_keys = set(domain.boundaries["vehicle"].cell_data.keys()) assert boundary_keys == {"C_p", "C_f", "normals"} def test_reserved_measure_key_cannot_be_a_target(self): with pytest.raises(ValueError, match="reserved"): - MeshToDomainMesh(cell_data_targets=[TARGET_QUADRATURE_MEASURE_KEY]) + MeshToDomainMesh(cell_data_targets=[EFFECTIVE_MEASURE_KEY]) @pytest.mark.parametrize("association", ["point_data", "cell_data"]) - def test_preexisting_reserved_measure_key_is_rejected(self, association): + def test_preexisting_effective_measures_are_preserved(self, association): mesh = _two_triangle_mesh_3d() - getattr(mesh, association)[TARGET_QUADRATURE_MEASURE_KEY] = torch.ones( + getattr(mesh, association)[EFFECTIVE_MEASURE_KEY] = torch.ones( mesh.n_points if association == "point_data" else mesh.n_cells ) - with pytest.raises(ValueError, match="already contains reserved key"): - MeshToDomainMesh(cell_data_targets=["C_p"])(mesh) + domain = MeshToDomainMesh(cell_data_targets=["C_p"])(mesh) + torch.testing.assert_close( + getattr(domain.boundaries["vehicle"], association)[EFFECTIVE_MEASURE_KEY], + getattr(mesh, association)[EFFECTIVE_MEASURE_KEY], + ) + torch.testing.assert_close( + domain.interior.point_data[EFFECTIVE_MEASURE_KEY], cell_measures(mesh) + ) def test_custom_boundary_name(self): mesh = _two_triangle_mesh_3d() diff --git a/test/datapipes/transforms/test_nested_fields.py b/test/datapipes/transforms/test_nested_fields.py index 473d132797..f3d8bf89eb 100644 --- a/test/datapipes/transforms/test_nested_fields.py +++ b/test/datapipes/transforms/test_nested_fields.py @@ -29,7 +29,6 @@ import physicsnemo.datapipes as dp from physicsnemo.datapipes.transforms.mesh import ( - TARGET_QUADRATURE_MEASURE_KEY, ComputeSurfaceNormals, DropMeshFields, MeshToDomainMesh, @@ -39,6 +38,7 @@ SetGlobalField, ) from physicsnemo.mesh import DomainMesh, Mesh +from physicsnemo.mesh.calculus.measure import EFFECTIVE_MEASURE_KEY def _surface_mesh() -> Mesh: @@ -280,7 +280,7 @@ def test_nested_target_moved_to_interior(self): domain = MeshToDomainMesh(cell_data_targets=["solution.pMeanTrim"])(mesh) assert _leaves(domain.interior.point_data) == { ("solution", "pMeanTrim"), - TARGET_QUADRATURE_MEASURE_KEY, + EFFECTIVE_MEASURE_KEY, } boundary = domain.boundaries["vehicle"] assert ("solution", "pMeanTrim") not in boundary.cell_data diff --git a/test/examples/test_globe_measures.py b/test/examples/test_globe_measures.py new file mode 100644 index 0000000000..e7df9b1d29 --- /dev/null +++ b/test/examples/test_globe_measures.py @@ -0,0 +1,62 @@ +# SPDX-FileCopyrightText: Copyright (c) 2023 - 2026 NVIDIA CORPORATION & AFFILIATES. +# SPDX-FileCopyrightText: All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Measure conservation when GLOBE resamples an already weighted mesh.""" + +import importlib.util +from pathlib import Path + +import pytest +import torch + +from physicsnemo.mesh import Mesh +from physicsnemo.mesh.calculus.measure import cell_measures, set_cell_measures + + +@pytest.mark.parametrize("geometry_only", [False, True]) +@pytest.mark.parametrize("explicit", [False, True]) +def test_subsample_preserves_represented_measure(geometry_only, explicit): + pytest.importorskip("pyvista") + path = ( + Path(__file__).parents[2] + / "examples/cfd/external_aerodynamics/globe/drivaer/dataset.py" + ) + spec = importlib.util.spec_from_file_location("globe_drivaer_dataset", path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + + points = torch.tensor([[0.0, 0, 0], [2, 0, 0], [0, 1, 0]]).repeat(8, 1) + mesh = Mesh( + points=points, + cells=torch.arange(24).reshape(8, 3), + cell_data={"id": torch.arange(8)}, + ) + if explicit: + set_cell_measures(mesh, torch.arange(1, 9, dtype=torch.float32)) + original = cell_measures(mesh).clone() + for n_cells in (4, 2): + torch.manual_seed(17) + indices = torch.randperm(mesh.n_cells)[:n_cells] + retained = cell_measures(mesh)[indices] + torch.manual_seed(17) + mesh = module.DrivAerMLDataSet.subsample_mesh( + mesh, n_cells, geometry_only=geometry_only + ) + torch.testing.assert_close(cell_measures(mesh).sum(), original.sum()) + torch.testing.assert_close( + cell_measures(mesh) / cell_measures(mesh).sum(), + retained / retained.sum(), + ) diff --git a/test/mesh/calculus/test_measure.py b/test/mesh/calculus/test_measure.py index d6b0d228bc..0fa9ba9d72 100644 --- a/test/mesh/calculus/test_measure.py +++ b/test/mesh/calculus/test_measure.py @@ -14,16 +14,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -"""Tests for ``physicsnemo.mesh.calculus.measure`` and weighted integration. - -Covers the module contract (ones fallback, multiplicative composition, -storage in ``cell_data`` under the reserved key), that ``integrate`` / -``integrate_flux`` / ``integrate_moment`` consume the effective measure -``cell_measures = cell_areas * measure_weights``, that a -Horvitz-Thompson-weighted cell subsample yields unbiased integrals, and -that weights survive slicing and rigid/scaling transforms with the -correct semantics. -""" +"""Effective cell measures, composable sampling corrections, and integration.""" import math @@ -33,10 +24,9 @@ from physicsnemo.mesh import Mesh from physicsnemo.mesh.calculus import integrate_moment from physicsnemo.mesh.calculus.measure import ( - MEASURE_WEIGHTS_KEY, - cell_measure_weights, + EFFECTIVE_MEASURE_KEY, cell_measures, - compose_measure_weights, + scale_measures, ) from physicsnemo.mesh.primitives.basic import two_triangles_2d @@ -67,53 +57,53 @@ def make_triangle_strip(n_cells: int, widths: torch.Tensor | None = None) -> Mes class TestSamplingWeights: def test_defaults_to_ones(self): mesh = two_triangles_2d.load() - w = cell_measure_weights(mesh) + w = cell_measures(mesh) assert w.shape == (mesh.n_cells,) - torch.testing.assert_close(w, torch.ones(mesh.n_cells)) + torch.testing.assert_close(w, mesh.cell_areas) ### The fallback must not materialize the reserved key. - assert MEASURE_WEIGHTS_KEY not in mesh.cell_data.keys() + assert EFFECTIVE_MEASURE_KEY not in mesh.cell_data.keys() ### And the effective measure is exactly the geometric one. torch.testing.assert_close(cell_measures(mesh), mesh.cell_areas) def test_compose_roundtrip_via_reserved_key(self): mesh = two_triangles_2d.load() - compose_measure_weights(mesh, torch.tensor([2.0, 3.0])) - assert MEASURE_WEIGHTS_KEY in mesh.cell_data.keys() - torch.testing.assert_close(cell_measure_weights(mesh), torch.tensor([2.0, 3.0])) + scale_measures(mesh, torch.tensor([2.0, 3.0])) + assert EFFECTIVE_MEASURE_KEY in mesh.cell_data.keys() + torch.testing.assert_close( + cell_measures(mesh), mesh.cell_areas * torch.tensor([2.0, 3.0]) + ) ### Stages compose multiplicatively. - compose_measure_weights(mesh, 10.0) + scale_measures(mesh, 10.0) torch.testing.assert_close( - cell_measure_weights(mesh), torch.tensor([20.0, 30.0]) + cell_measures(mesh), mesh.cell_areas * torch.tensor([20.0, 30.0]) ) def test_storage_rejects_wrong_shape(self): ### cell_data's batch dimension rejects a wrong leading dimension. mesh = two_triangles_2d.load() with pytest.raises(RuntimeError): - mesh.cell_data[MEASURE_WEIGHTS_KEY] = torch.ones(mesh.n_cells + 1) + mesh.cell_data[EFFECTIVE_MEASURE_KEY] = torch.ones(mesh.n_cells + 1) def test_reserved_field_rejects_trailing_singleton_dimension(self): ### TensorDict accepts vector-valued cell data, so the reserved ### scalar field must enforce its own exact shape. mesh = two_triangles_2d.load() - mesh.cell_data[MEASURE_WEIGHTS_KEY] = torch.ones(mesh.n_cells, 1) + mesh.cell_data[EFFECTIVE_MEASURE_KEY] = torch.ones(mesh.n_cells, 1) - with pytest.raises(ValueError, match="one scalar per cell"): - cell_measure_weights(mesh) with pytest.raises(ValueError, match="one scalar per cell"): cell_measures(mesh) def test_compose_rejects_non_scalar_broadcast_shape(self): mesh = two_triangles_2d.load() with pytest.raises(ValueError, match="scalar or have shape"): - compose_measure_weights(mesh, torch.ones(mesh.n_cells, 1)) + scale_measures(mesh, torch.ones(mesh.n_cells, 1)) def test_weights_survive_slice_cells(self): mesh = make_triangle_strip(6) - compose_measure_weights(mesh, torch.arange(1.0, 7.0)) + scale_measures(mesh, torch.arange(1.0, 7.0)) sliced = mesh.slice_cells(torch.tensor([1, 4])) torch.testing.assert_close( - cell_measure_weights(sliced), torch.tensor([2.0, 5.0]) + cell_measures(sliced), sliced.cell_areas * torch.tensor([2.0, 5.0]) ) @@ -122,14 +112,14 @@ def test_integrate_cell_data_uses_effective_measure(self): mesh = make_triangle_strip(4) mesh.cell_data["f"] = torch.tensor([1.0, 2.0, 3.0, 4.0]) unweighted = mesh.integrate("f") - compose_measure_weights(mesh, torch.full((4,), 2.5)) + scale_measures(mesh, torch.full((4,), 2.5)) torch.testing.assert_close(mesh.integrate("f"), unweighted * 2.5) def test_integrate_point_data_uses_effective_measure(self): mesh = make_triangle_strip(3) mesh.point_data["T"] = torch.randn(mesh.n_points) unweighted = mesh.integrate("T", data_source="points") - compose_measure_weights(mesh, torch.full((3,), 4.0)) + scale_measures(mesh, torch.full((3,), 4.0)) torch.testing.assert_close( mesh.integrate("T", data_source="points"), unweighted * 4.0 ) @@ -138,7 +128,7 @@ def test_integrate_flux_uses_effective_measure(self): mesh = make_triangle_strip(3) # planar, normals +/- z mesh.cell_data["v"] = torch.randn(3, 3) unweighted = mesh.integrate_flux("v") - compose_measure_weights(mesh, torch.full((3,), 3.0)) + scale_measures(mesh, torch.full((3,), 3.0)) torch.testing.assert_close(mesh.integrate_flux("v"), unweighted * 3.0) def test_integrate_moment_uses_effective_measure(self): @@ -146,7 +136,7 @@ def test_integrate_moment_uses_effective_measure(self): left = torch.randn(4, 2) right = torch.randn(4, 3) unweighted = integrate_moment(mesh, left, right) - compose_measure_weights(mesh, torch.full((4,), 2.0)) + scale_measures(mesh, torch.full((4,), 2.0)) torch.testing.assert_close( integrate_moment(mesh, left, right), unweighted * 2.0 ) @@ -171,7 +161,7 @@ def test_ht_subsample_integral_unbiased_over_all_starts(self): for start in range(n): idx = torch.arange(start, start + k) % n sub = mesh.slice_cells(idx) - compose_measure_weights(sub, torch.full((k,), n / k)) + scale_measures(sub, torch.full((k,), n / k)) estimates.append(sub.integrate("f").to(torch.float64)) torch.testing.assert_close( torch.stack(estimates).mean(), full, rtol=1e-5, atol=1e-6 @@ -186,7 +176,7 @@ def test_weights_invariant_under_rigid_and_scaling_transforms(self): n = 5 mesh = make_triangle_strip(n, widths=torch.rand(n) + 0.5) mesh.cell_data["f"] = torch.ones(n) - compose_measure_weights(mesh, torch.full((n,), 2.0)) + scale_measures(mesh, torch.full((n,), 2.0)) base_integral = mesh.integrate("f") moved = ( @@ -195,7 +185,7 @@ def test_weights_invariant_under_rigid_and_scaling_transforms(self): .scale(1.0 / 5.0, transform_cell_data=True) ) - torch.testing.assert_close(cell_measure_weights(moved), torch.full((n,), 2.0)) + torch.testing.assert_close(cell_measures(moved), moved.cell_areas * 2.0) torch.testing.assert_close( moved.integrate("f"), base_integral / 25.0, rtol=1e-5, atol=1e-7 ) diff --git a/test/mesh/calculus/test_point_measures.py b/test/mesh/calculus/test_point_measures.py new file mode 100644 index 0000000000..5372a3ff04 --- /dev/null +++ b/test/mesh/calculus/test_point_measures.py @@ -0,0 +1,324 @@ +# SPDX-FileCopyrightText: Copyright (c) 2023 - 2026 NVIDIA CORPORATION & AFFILIATES. +# SPDX-FileCopyrightText: All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Contracts for explicit measures across point and cell representations.""" + +import pytest +import torch + +from physicsnemo.datapipes.readers.mesh import _subsample_mesh_points +from physicsnemo.datapipes.transforms.mesh import ( + MeshToDomainMesh, + ScaleMesh, + SubsampleMesh, +) +from physicsnemo.mesh import Mesh +from physicsnemo.mesh.calculus.measure import ( + EFFECTIVE_MEASURE_KEY, + POINT_MEASURE_DIMENSION_KEY, + cell_measures, + lumped_point_measures, + point_measures, + scale_measures, + set_cell_measures, + set_point_measures, +) +from physicsnemo.mesh.primitives.basic import two_triangles_2d + + +def simplex(dimension, *, spatial_dimension=3): + points = torch.cat( + [torch.zeros(1, spatial_dimension), torch.eye(spatial_dimension)[:dimension]] + ) + return Mesh(points=points, cells=torch.arange(dimension + 1).unsqueeze(0)) + + +@pytest.mark.parametrize("connected", [False, True]) +def test_point_quadrature_never_infers_a_measure(connected): + source = two_triangles_2d.load() + mesh = source if connected else Mesh(points=source.points) + with pytest.raises(KeyError, match="explicit effective measures"): + point_measures(mesh) + with pytest.raises(KeyError, match="explicit effective measures"): + mesh.integrate_samples(torch.ones(mesh.n_points)) + + values = torch.arange(mesh.n_points, dtype=mesh.points.dtype) + 1 + weights = torch.arange(mesh.n_points, dtype=mesh.points.dtype) + 2 + set_point_measures(mesh, weights, dimension=2) + torch.testing.assert_close(mesh.integrate_samples(values), (values * weights).sum()) + assert EFFECTIVE_MEASURE_KEY not in mesh.cell_data + + +def test_nodal_and_sample_integration_are_explicitly_distinct(): + mesh = two_triangles_2d.load() + scale_measures(mesh, torch.tensor([2.0, 3.0])) + values = torch.arange(mesh.n_points, dtype=mesh.points.dtype) + p1 = mesh.integrate(values, data_source="points") + set_point_measures(mesh, torch.ones(mesh.n_points), dimension=0) + torch.testing.assert_close(mesh.integrate_samples(values), values.sum()) + torch.testing.assert_close(mesh.integrate(values, data_source="points"), p1) + set_point_measures(mesh, lumped_point_measures(mesh), dimension=2) + torch.testing.assert_close(mesh.integrate_samples(values), p1) + torch.testing.assert_close(point_measures(mesh).sum(), cell_measures(mesh).sum()) + + +@pytest.mark.parametrize("dimension", [1, 2, 3]) +@pytest.mark.parametrize("factor", [-3.0, 0.0, 2.0]) +@pytest.mark.parametrize("transform_fields", [False, True]) +def test_centroid_conversion_commutes_with_uniform_scaling( + dimension, factor, transform_fields +): + source = simplex(dimension) + source.cell_data["target"] = torch.tensor([7.0]) + scale_measures(source, 2.5) + convert = MeshToDomainMesh(cell_data_targets=["target"]) + domain = convert(source) + moved = ScaleMesh( + factor, + transform_point_data=transform_fields, + transform_cell_data=transform_fields, + ).apply_to_domain(domain) + expected = cell_measures(source) * abs(factor) ** dimension + torch.testing.assert_close(point_measures(moved.interior), expected) + torch.testing.assert_close(cell_measures(moved.boundaries["vehicle"]), expected) + torch.testing.assert_close( + point_measures(convert(source.scale(factor)).interior), expected + ) + torch.testing.assert_close(point_measures(domain.interior), cell_measures(source)) + assert int(moved.interior.global_data[POINT_MEASURE_DIMENSION_KEY]) == dimension + + +def test_cell_measures_follow_anisotropic_deformation(): + mesh = simplex(2) + scale_measures(mesh, 3.0) + moved = mesh.scale([2.0, 3.0, 4.0]) + torch.testing.assert_close(cell_measures(moved), moved.cell_areas * 3.0) + displaced = mesh.with_points(mesh.points * torch.tensor([2.0, 3.0, 4.0])) + torch.testing.assert_close(cell_measures(displaced), cell_measures(moved)) + + +def test_point_support_is_required_for_anisotropic_surface_scaling(): + domain = MeshToDomainMesh()(simplex(2)) + before = point_measures(domain.interior).clone() + with pytest.raises(ValueError, match="support geometry"): + domain.scale([2.0, 3.0, 4.0]) + with pytest.raises(ValueError, match="preservation policy"): + domain.interior.with_points(domain.interior.points * 2) + fixed = domain.interior.with_points( + domain.interior.points * 2, preserve_measures=True + ) + torch.testing.assert_close(point_measures(fixed), before) + torch.testing.assert_close(point_measures(domain.interior), before) + + +def test_counting_and_volume_point_measures_have_known_affine_scaling(): + cloud = Mesh(points=torch.randn(5, 3)) + set_point_measures(cloud, torch.ones(5), dimension=0) + torch.testing.assert_close( + point_measures(cloud.scale([2.0, 3.0, 4.0])), torch.ones(5) + ) + torch.testing.assert_close( + point_measures(cloud.with_points(cloud.points + 1)), torch.ones(5) + ) + set_point_measures(cloud, torch.full((5,), 2.0), dimension=3) + torch.testing.assert_close( + point_measures(cloud.scale([2.0, 3.0, 4.0])), torch.full((5,), 48.0) + ) + + +def test_rigid_transforms_and_dtype_preserve_point_measures(): + cloud = MeshToDomainMesh()(simplex(2)).interior + moved = ( + cloud.rotate(0.7, axis="x").translate([1.0, 2.0, 3.0]).to(dtype=torch.float64) + ) + torch.testing.assert_close(point_measures(moved), point_measures(cloud).double()) + torch.testing.assert_close( + point_measures(moved.scale(2.0)), point_measures(cloud).double() * 4 + ) + + +@pytest.mark.parametrize("filter", ["linear", "loop", "butterfly"]) +def test_subdivision_distributes_complete_cell_measures(filter): + mesh = two_triangles_2d.load() + scale_measures(mesh, 3.0) + refined = mesh.subdivide(levels=2, filter=filter) + torch.testing.assert_close(cell_measures(refined), refined.cell_areas * 3.0) + if filter == "linear": + torch.testing.assert_close( + cell_measures(refined).sum(), cell_measures(mesh).sum() + ) + set_point_measures(mesh, torch.ones(mesh.n_points), dimension=0) + with pytest.raises(ValueError, match="cannot be interpolated"): + mesh.subdivide(filter=filter) + + +@pytest.mark.parametrize("reader", [False, True]) +def test_point_sampling_reweights_only_explicit_quadrature(reader): + cloud = Mesh(points=torch.randn(10, 3)) + set_point_measures(cloud, torch.ones(10), dimension=0) + if reader: + sampled = _subsample_mesh_points(cloud, 4, torch.Generator().manual_seed(1)) + else: + sampled = SubsampleMesh(n_points=4)(cloud) + torch.testing.assert_close(point_measures(sampled), torch.full((4,), 2.5)) + torch.testing.assert_close(point_measures(cloud), torch.ones(10)) + ordinary = SubsampleMesh(n_points=4)(Mesh(points=cloud.points)) + with pytest.raises(KeyError): + point_measures(ordinary) + + +def test_measure_storage_survives_serialization_slice_and_merge(tmp_path): + cloud = Mesh(points=torch.randn(5, 3), point_data={"value": torch.arange(5.0)}) + set_point_measures(cloud, torch.arange(5.0) + 1, dimension=2) + cloud.save(tmp_path / "cloud.pmsh") + loaded = Mesh.load(tmp_path / "cloud.pmsh") + selected = loaded.slice_points([1, 3]) + torch.testing.assert_close(point_measures(selected), torch.tensor([2.0, 4.0])) + merged = Mesh.merge([selected, selected]) + torch.testing.assert_close( + point_measures(merged.scale(2.0)), torch.tensor([8.0, 16.0, 8.0, 16.0]) + ) + assert merged.global_data[POINT_MEASURE_DIMENSION_KEY].shape == () + incompatible = selected.clone() + set_point_measures(incompatible, point_measures(incompatible), dimension=1) + with pytest.raises(ValueError, match="different measure dimensions"): + Mesh.merge([selected, incompatible]) + + +def test_centroid_conversion_and_vertex_conversion_preserve_existing_measures(): + source = simplex(2) + set_cell_measures(source, torch.tensor([5.0])) + set_point_measures(source, torch.tensor([1.0, 2.0, 3.0]), dimension=2) + vertices = MeshToDomainMesh(interior_points="vertices")(source) + torch.testing.assert_close( + point_measures(vertices.interior), point_measures(source) + ) + centroids = MeshToDomainMesh()(source) + torch.testing.assert_close(point_measures(centroids.interior), torch.tensor([5.0])) + torch.testing.assert_close( + point_measures(centroids.boundaries["vehicle"]), point_measures(source) + ) + + +@pytest.mark.parametrize("explicit", [False, True]) +@pytest.mark.parametrize("dual", [False, True]) +def test_centroid_mesh_conversions_transfer_measures_without_mutating_source( + explicit, dual +): + """Centroid quadrature keeps the source dimension, including on dual graphs.""" + source = two_triangles_2d.load() + source.cell_data["value"] = torch.tensor([2.0, 5.0]) + if explicit: + scale_measures(source, torch.tensor([3.0, 7.0])) + # Source vertex quadrature is independent of the new centroid quadrature. + set_point_measures(source, torch.ones(source.n_points), dimension=0) + expected = cell_measures(source).clone() + converted = ( + source.to_dual_graph() + if dual + else source.to_point_cloud(point_source="cell_centroids") + ) + torch.testing.assert_close(point_measures(converted), expected) + torch.testing.assert_close( + converted.integrate_samples("value"), source.integrate("value") + ) + torch.testing.assert_close(point_measures(converted.scale(2.0)), expected * 4) + assert int(converted.global_data[POINT_MEASURE_DIMENSION_KEY]) == 2 + assert int(source.global_data[POINT_MEASURE_DIMENSION_KEY]) == 0 + assert (EFFECTIVE_MEASURE_KEY in source.cell_data) == explicit + if dual: + torch.testing.assert_close(cell_measures(converted), converted.cell_areas) + scale_measures(converted, 2.0, association="points") + torch.testing.assert_close(cell_measures(source), expected) + + +def test_generic_field_interpolation_does_not_interpolate_measures(): + mesh = simplex(2) + set_cell_measures(mesh, torch.tensor([7.0])) + mesh.cell_data["f"] = torch.tensor([3.0]) + converted = mesh.cell_data_to_point_data() + assert EFFECTIVE_MEASURE_KEY not in converted.point_data + torch.testing.assert_close(converted.point_data["f"], torch.full((3,), 3.0)) + set_point_measures(mesh, torch.ones(3), dimension=2) + mesh.point_data["g"] = torch.tensor([1.0, 2.0, 3.0]) + converted = mesh.point_data_to_cell_data() + torch.testing.assert_close(cell_measures(converted), torch.tensor([7.0])) + torch.testing.assert_close(converted.cell_data["g"], torch.tensor([2.0])) + + +def test_measure_validation_rejects_ambiguous_shapes_and_dimensions(): + mesh = simplex(2) + with pytest.raises(ValueError, match="one scalar"): + set_cell_measures(mesh, torch.ones(1, 1)) + with pytest.raises(ValueError, match="one scalar"): + set_point_measures(mesh, torch.ones(3, 1), dimension=2) + with pytest.raises(ValueError, match="dimension"): + set_point_measures(mesh, torch.ones(3), dimension=4) + set_point_measures(mesh, torch.ones(3), dimension=2) + before = point_measures(mesh).clone() + with pytest.raises(ValueError, match="scalar or have shape"): + scale_measures(mesh, torch.ones(3, 1), association="points") + torch.testing.assert_close(point_measures(mesh), before) + + +def test_differentiable_sample_quadrature(): + cloud = Mesh(points=torch.zeros(3, 2)) + values = torch.tensor([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]], requires_grad=True) + measures = torch.tensor([2.0, 3.0, 4.0], requires_grad=True) + set_point_measures(cloud, measures, dimension=2) + cloud.integrate_samples(values).sum().backward() + torch.testing.assert_close( + values.grad, measures.detach()[:, None].expand_as(values) + ) + torch.testing.assert_close(measures.grad, values.detach().sum(-1)) + + +def test_empty_sample_integral_and_nan_omission(): + cloud = Mesh(points=torch.empty(0, 3)) + set_point_measures(cloud, torch.empty(0), dimension=2) + torch.testing.assert_close( + cloud.integrate_samples(torch.empty(0, 2)), torch.zeros(2) + ) + cloud = Mesh(points=torch.zeros(3, 2)) + set_point_measures(cloud, torch.tensor([1.0, 2.0, 3.0]), dimension=2) + values = torch.tensor([1.0, float("nan"), 4.0]) + torch.testing.assert_close(cloud.integrate_samples(values), torch.tensor(13.0)) + assert torch.isnan(cloud.integrate_samples(values, nan_policy="propagate")) + + +def test_legacy_multiplier_cannot_silently_become_geometric_fallback(): + mesh = simplex(2) + mesh.cell_data["_measure_weights"] = torch.tensor([3.0]) + with pytest.raises(ValueError, match="must be converted"): + cell_measures(mesh) + + +def test_remeshing_requires_a_measure_transfer_rule(): + mesh = simplex(2) + set_cell_measures(mesh, torch.tensor([1.0])) + with pytest.raises(ValueError, match="conservative measure transfer"): + mesh.remesh(n_clusters=3) + + +def test_partition_accumulates_effective_measures_without_overriding_geometry(): + from physicsnemo.mesh.remeshing import partition_cells + + mesh = two_triangles_2d.load() + geometric = mesh.cell_areas.clone() + set_cell_measures(mesh, torch.tensor([2.0, 5.0])) + partition = partition_cells(mesh, seeds=mesh.cell_centroids) + torch.testing.assert_close(partition.cluster_areas, torch.tensor([2.0, 5.0])) + torch.testing.assert_close(mesh.cell_areas, geometric) diff --git a/test/mesh/io/io_zarr/test_reader_integration.py b/test/mesh/io/io_zarr/test_reader_integration.py index 3000148b4e..f4f0861edd 100644 --- a/test/mesh/io/io_zarr/test_reader_integration.py +++ b/test/mesh/io/io_zarr/test_reader_integration.py @@ -20,6 +20,7 @@ from pathlib import Path +import pytest import torch from conftest import assert_meshes_equal, make_domain_mesh, make_mesh @@ -28,6 +29,8 @@ MeshReader, _subsample_mesh, ) +from physicsnemo.mesh import DomainMesh, Mesh +from physicsnemo.mesh.calculus.measure import point_measures, set_point_measures from physicsnemo.mesh.io import from_zarr, to_zarr @@ -109,3 +112,31 @@ def test_mixed_directory_discovery(tmp_path): metas = {Path(reader[i][1]["source_path"]).name: reader[i][0] for i in (0, 1)} assert_meshes_equal(m1, metas["a.mesh.zarr"]) assert_meshes_equal(m2, metas["b.pmsh"]) + + +@pytest.mark.parametrize("domain", [False, True]) +@pytest.mark.parametrize("n_keep", [4, 20]) +def test_point_quadrature_subsampling_matches_between_formats(tmp_path, domain, n_keep): + """Partial reads apply the same correction as eager reads, exactly once.""" + cloud = Mesh( + points=torch.arange(30, dtype=torch.float32).reshape(10, 3), + point_data={"index": torch.arange(10)}, + ) + set_point_measures(cloud, torch.arange(10, dtype=torch.float32) + 1, dimension=3) + sample = DomainMesh(interior=cloud) if domain else cloud + sample.save(tmp_path / "sample.pmsh") + to_zarr(sample, tmp_path / "sample.zarr", chunk_rows=2) + reader_type = DomainMeshReader if domain else MeshReader + loaded = [] + for pattern in ("*.pmsh", "*.zarr"): + reader = reader_type( + tmp_path, pattern=pattern, subsample_n_points=n_keep, pin_memory=False + ) + reader.set_generator(torch.Generator().manual_seed(5)) + sampled, _ = reader[0] + points = sampled.interior if domain else sampled + expected = (points.point_data["index"] + 1).float() * (10 / min(n_keep, 10)) + torch.testing.assert_close(point_measures(points), expected) + torch.testing.assert_close(point_measures(points.scale(2.0)), expected * 8) + loaded.append(points) + assert_meshes_equal(*loaded) diff --git a/test/models/globe/test_measure.py b/test/models/globe/test_measure.py index 95758cca38..e478b101ca 100644 --- a/test/models/globe/test_measure.py +++ b/test/models/globe/test_measure.py @@ -17,7 +17,7 @@ """GLOBE consumption of the effective cell measure. GLOBE weights its boundary integrals by the effective cell measure -``cell_areas * measure_weights`` (see +``cell_measures(mesh)`` (see :mod:`physicsnemo.mesh.calculus.measure`), and compounds one such measure-weighted sum per integral stage (``n_communication_hyperlayers + 1`` in total), so an incorrect measure is @@ -35,7 +35,7 @@ import torch from physicsnemo.experimental.models.globe.model import GLOBE -from physicsnemo.mesh.calculus.measure import compose_measure_weights +from physicsnemo.mesh.calculus.measure import scale_measures from physicsnemo.mesh.primitives.procedural import lumpy_sphere SEED = 7 @@ -92,7 +92,7 @@ def test_weights_equivalent_to_scaled_areas(): for name in ("vehicle", "floor"): mesh_w = kwargs_weighted["boundary_meshes"][name] w = torch.rand(mesh_w.n_cells, generator=gen) + 0.5 - compose_measure_weights(mesh_w, w) + scale_measures(mesh_w, w) mesh_s = kwargs_scaled["boundary_meshes"][name] mesh_s._cache["cell", "areas"] = mesh_s.cell_areas * w @@ -112,7 +112,7 @@ def test_unit_weights_match_no_weights(): kwargs_ones = _make_inputs(device) for name in ("vehicle", "floor"): mesh = kwargs_ones["boundary_meshes"][name] - compose_measure_weights(mesh, torch.ones(mesh.n_cells)) + scale_measures(mesh, torch.ones(mesh.n_cells)) out_plain = _forward(model, kwargs_plain) out_ones = _forward(model, kwargs_ones) @@ -134,7 +134,7 @@ def test_weights_change_output(): kwargs_weighted = _make_inputs(device) for name in ("vehicle", "floor"): mesh = kwargs_weighted["boundary_meshes"][name] - compose_measure_weights(mesh, torch.full((mesh.n_cells,), 3.0)) + scale_measures(mesh, torch.full((mesh.n_cells,), 3.0)) out_plain = _forward(model, kwargs_plain) out_weighted = _forward(model, kwargs_weighted) From 241a3e9f2259a1f14973d0f52f27e437ff6e395e Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Sat, 19 Sep 2026 16:40:58 -0400 Subject: [PATCH 5/6] Keep GLOBE dataset measure tests with the GLOBE suite --- .../globe/test_dataset_measure.py} | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename test/{examples/test_globe_measures.py => models/globe/test_dataset_measure.py} (98%) diff --git a/test/examples/test_globe_measures.py b/test/models/globe/test_dataset_measure.py similarity index 98% rename from test/examples/test_globe_measures.py rename to test/models/globe/test_dataset_measure.py index e7df9b1d29..517255420c 100644 --- a/test/examples/test_globe_measures.py +++ b/test/models/globe/test_dataset_measure.py @@ -31,7 +31,7 @@ def test_subsample_preserves_represented_measure(geometry_only, explicit): pytest.importorskip("pyvista") path = ( - Path(__file__).parents[2] + Path(__file__).parents[3] / "examples/cfd/external_aerodynamics/globe/drivaer/dataset.py" ) spec = importlib.util.spec_from_file_location("globe_drivaer_dataset", path) From e2ff58109108d67da4bb26eb8bb35a47b6cb045e Mon Sep 17 00:00:00 2001 From: Peter Sharpe Date: Sat, 19 Sep 2026 19:08:51 -0400 Subject: [PATCH 6/6] Tighten mesh calculus documentation and clarify released migration --- CHANGELOG.md | 17 +++++ docs/api/mesh/calculus.rst | 131 +++++++------------------------------ 2 files changed, 42 insertions(+), 106 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c6ec81fd71..9da92bdcbb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,6 +44,23 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 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 diff --git a/docs/api/mesh/calculus.rst b/docs/api/mesh/calculus.rst index 0a48170c43..0a4cfff3bc 100644 --- a/docs/api/mesh/calculus.rst +++ b/docs/api/mesh/calculus.rst @@ -53,120 +53,39 @@ Key Operators Effective measures and integration ---------------------------------- -The reserved ``_effective_measure`` field stores the complete measure associated -with each cell or point sample. It has shape ``(n_cells,)`` in ``cell_data`` or -``(n_points,)`` in ``point_data``. It is never a dimensionless correction that -must still be multiplied by a geometric area or volume. +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. -``cell_measures(mesh)`` returns the explicit cell measures, or the geometric -simplex measures when the field is absent. ``point_measures(mesh)`` requires -explicit point measures: it never changes its interpretation based on whether -``mesh.cells`` is empty. An ordinary sum is counting measure; it can also be -represented explicitly by installing ones with ``dimension=0``. - -.. code:: python - - from physicsnemo.mesh.calculus import ( - cell_measures, - point_measures, - scale_measures, - set_point_measures, - ) - - # A sampling stage retains k of N cells. Prior corrections are preserved. - sampled = mesh.slice_cells(indices) - scale_measures(sampled, mesh.n_cells / sampled.n_cells) +* ``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. - # Transfer complete measures when cells become independent point samples. - queries = Mesh(points=sampled.cell_centroids) - set_point_measures( - queries, cell_measures(sampled), dimension=sampled.n_manifold_dims - ) - integral = queries.integrate_samples(predictions) +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. -There are two distinct integration operations: - -* ``mesh.integrate(field, data_source="cells")`` integrates piecewise-constant - cell values over the cells. ``data_source="points"`` integrates a - piecewise-linear vertex field over those same cells. Both use effective - **cell** measures. Vertex fields keep the existing rule that a NaN vertex - invalidates its incident cell's contribution when ``nan_policy="omit"``. -* ``mesh.integrate_samples(field)`` sums independent point samples times explicit - **point** measures, with no dependence on connectivity. NaN omission applies - independently to each sample contribution. Missing point measures raise an - error; there is no implicit counting or geometric fallback. - -``lumped_point_measures(mesh)`` explicitly constructs vertex quadrature by -sharing each cell's effective measure equally among its vertices. For finite -nodal fields it reproduces piecewise-linear integration. Calling it does not -modify the mesh or change either integration rule. - -Measure lifecycle -~~~~~~~~~~~~~~~~~ - -All storage and reweighting helpers live in ``physicsnemo.mesh.calculus.measure``. -``set_cell_measures`` and ``set_point_measures`` assign complete measures; -``scale_measures`` multiplies existing measures by a scalar or per-entity factor. -Sampling uses the latter with inverse inclusion probabilities. Raw slicing is a -restriction and does not apply a sampling correction. Serialization and device -transfers preserve the fields. - -``mesh.to_point_cloud(point_source="cell_centroids")`` and ``mesh.to_dual_graph`` -transfer each cell's complete measure to its centroid, including the represented -dimension. Dual-graph edges retain their own geometric length measure; the -original surface or volume measure is associated with the graph's points. - -Point measures carry a scalar ``_point_measure_dimension`` in their mesh's -``global_data``. ``set_point_measures`` writes this metadata: 0 for counting, 1 -for length, 2 for area, and 3 for volume. It describes the represented measure, -not the point cloud's topological dimension. Merging point quadrature requires -matching measure dimensions and preserves this scalar metadata. - -Rigid transformations preserve measures. Uniform scaling by ``s`` multiplies -measures of dimension ``d`` by ``abs(s)**d``. Cell geometry changes preserve the -ratio of represented to geometric measure. Subdivision transfers that ratio to -children; linear subdivision therefore conserves each parent's total measure. -A nonzero measure on a geometrically degenerate cell needs explicit replacement -measures when its geometry changes. - -For points, full-dimensional measures also support general square linear maps -through their absolute determinant. Anisotropic transformations of embedded -surface/curve quadrature require support geometry that a point cloud does not -contain, and are rejected. Transform the source cells before creating those -samples, or explicitly retain reference measures with -``mesh.with_points(new_points, preserve_measures=True)``. This policy also serves -coordinate normalization where measures deliberately remain in reference units. - -Generic field interpolation excludes effective measures. Subdivision of explicit -point quadrature and remeshing of explicit cell/point quadrature require an -explicit conservative transfer or replacement measures; ordinary field -interpolation cannot provide one. - -This pre-release API replaces ``_measure_weights`` and the datapipes-specific -target quadrature key. Files using the old cell multiplier must be regenerated -or converted to complete measures before use. Existing geometric areas remain -geometric: producers must not override area caches to store represented areas. - -For a legacy mesh, convert the stored multiplier once and remove its old key: +Converting cells to centroid samples transfers their measures automatically: .. code:: python - from physicsnemo.mesh.calculus import set_cell_measures, set_point_measures - - if "_measure_weights" in mesh.cell_data: - weights = mesh.cell_data.pop("_measure_weights") - set_cell_measures(mesh, mesh.cell_areas * weights) + queries = mesh.to_point_cloud(point_source="cell_centroids") + values = queries.points[:, 0] # Integrate f(x, ...) = x. + integral = queries.integrate_samples(values) - # Legacy centroid query measures already contain the geometric contribution. - # Use the dimension of the source cells (2 here for a surface), not the - # zero-dimensional topology of the query cloud. - if "_target_quadrature_measure" in mesh.point_data: - measures = mesh.point_data.pop("_target_quadrature_measure") - set_point_measures(mesh, measures, dimension=2) +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. -Update field mappings to ``cell_data._effective_measure`` or -``point_data._effective_measure`` as appropriate. Do not multiply these complete -measures by geometric areas a second time. +See :doc:`transformations` for supported geometric changes and explicit +preservation of reference measures. API Reference -------------