From ac46e469948f76f7b161356e2f0552a522338ada Mon Sep 17 00:00:00 2001 From: Emanuel Pituch <32016786+epituch@users.noreply.github.com> Date: Thu, 16 Apr 2026 20:29:47 -0700 Subject: [PATCH 1/3] Section 4 WIP --- pyproject.toml | 2 +- src/nmesh/backend.py | 18 +- src/nmesh/io/meshio_support.py | 7 +- src/nmesh/mesher/__init__.py | 1 + src/nmesh/mesher/meshing_parameters.py | 4 +- src/nmesh/mesher/periodic.py | 81 +++++++ src/nmesh/mesher/relaxation/__init__.py | 14 ++ src/nmesh/mesher/relaxation/_constants.py | 18 ++ src/nmesh/mesher/relaxation/_types.py | 14 ++ src/nmesh/mesher/relaxation/density.py | 131 ++++++++++ src/nmesh/mesher/relaxation/engine.py | 281 ++++++++++++++++++++++ src/nmesh/mesher/relaxation/geometry.py | 166 +++++++++++++ src/nmesh/mesher/relaxation/seeding.py | 225 +++++++++++++++++ src/nmesh/mesher/relaxation/topology.py | 193 +++++++++++++++ src/nmesh/utils/constants.py | 4 + tests/nmesh/test_meshing_engine.py | 46 ++++ tests/nmesh_test.py | 70 +++++- 17 files changed, 1257 insertions(+), 18 deletions(-) create mode 100644 src/nmesh/mesher/periodic.py create mode 100644 src/nmesh/mesher/relaxation/__init__.py create mode 100644 src/nmesh/mesher/relaxation/_constants.py create mode 100644 src/nmesh/mesher/relaxation/_types.py create mode 100644 src/nmesh/mesher/relaxation/density.py create mode 100644 src/nmesh/mesher/relaxation/engine.py create mode 100644 src/nmesh/mesher/relaxation/geometry.py create mode 100644 src/nmesh/mesher/relaxation/seeding.py create mode 100644 src/nmesh/mesher/relaxation/topology.py create mode 100644 tests/nmesh/test_meshing_engine.py diff --git a/pyproject.toml b/pyproject.toml index f67e694..348a1b8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,7 @@ requires-python = ">=3.10" classifiers = [ "Programming Language :: Python :: 3", "Operating System :: OS Independent" ] -dependencies = ["pint", "numpy", "tabulate", "meshio", "h5py"] +dependencies = ["pint", "numpy", "scipy", "tabulate", "meshio", "h5py"] [project.optional-dependencies] test = ["pytest", "pytest-cov", "pytest-watch"] diff --git a/src/nmesh/backend.py b/src/nmesh/backend.py index f24518d..e608bab 100644 --- a/src/nmesh/backend.py +++ b/src/nmesh/backend.py @@ -167,7 +167,23 @@ def mesh_bodies_raw( cache, hints, ): - return RawMesh(dim=len(bb_min)) + from .mesher.relaxation import mesh_bodies_raw as python_mesh_bodies_raw + + return python_mesh_bodies_raw( + driver, + mesher, + bb_min, + bb_max, + mesh_ext, + objects, + a0, + density, + fixed, + mobile, + simply, + periodic, + hints, + ) def mesh_from_points_and_simplices( self, diff --git a/src/nmesh/io/meshio_support.py b/src/nmesh/io/meshio_support.py index 33cffdb..a21bfa4 100644 --- a/src/nmesh/io/meshio_support.py +++ b/src/nmesh/io/meshio_support.py @@ -7,12 +7,7 @@ from ..backend import RawMesh from .legacy_nmesh_hdf5 import load_raw_mesh_from_legacy_nmesh_hdf5 - -try: - from meshio._exceptions import ReadError as MeshioReadError -except ImportError: - # Fallback for older meshio versions - MeshioReadError = Exception +from meshio._exceptions import ReadError as MeshioReadError _CELL_TYPE_BY_DIM = { diff --git a/src/nmesh/mesher/__init__.py b/src/nmesh/mesher/__init__.py index b1691f0..579ddf5 100644 --- a/src/nmesh/mesher/__init__.py +++ b/src/nmesh/mesher/__init__.py @@ -1,2 +1,3 @@ from .meshing_parameters import * from .driver import * +from .relaxation import * diff --git a/src/nmesh/mesher/meshing_parameters.py b/src/nmesh/mesher/meshing_parameters.py index 4451fa9..d366013 100644 --- a/src/nmesh/mesher/meshing_parameters.py +++ b/src/nmesh/mesher/meshing_parameters.py @@ -7,7 +7,7 @@ from mock_features import MockFeatures -from ..utils.constants import MIN_DIVISION_MAGNITUDE +from ..utils.constants import BOUNDARY_FUZZ, MIN_DIVISION_MAGNITUDE log = logging.getLogger(__name__) @@ -248,7 +248,7 @@ def _setup_defaults(self): # Volume determination "nr_probes_for_determining_volume": 100000, # Boundary condition parameters - "boundary_condition_acceptable_fuzz": 1e-6, + "boundary_condition_acceptable_fuzz": BOUNDARY_FUZZ, "boundary_condition_max_nr_correction_steps": 200, "boundary_condition_debuglevel": 0, # Relaxation parameters diff --git a/src/nmesh/mesher/periodic.py b/src/nmesh/mesher/periodic.py new file mode 100644 index 0000000..a4ecc04 --- /dev/null +++ b/src/nmesh/mesher/periodic.py @@ -0,0 +1,81 @@ +from __future__ import annotations + +"""Helpers for building periodic-equivalence groups from mesh coordinates.""" + +from collections import defaultdict +from collections.abc import Iterable + +import numpy as np + + +def _point_key(values: np.ndarray, decimals: int = 8) -> tuple[float, ...]: + """Return a rounded tuple key suitable for coordinate-based lookups.""" + + return tuple(np.round(values.astype(float, copy=False), decimals=decimals).tolist()) + + +def merge_periodic_groups(groups: Iterable[Iterable[int]]) -> list[list[int]]: + """Merge overlapping periodic point groups into disjoint sorted groups.""" + + merged: list[set[int]] = [] + for group in groups: + candidate = {int(index) for index in group} + if len(candidate) < 2: + continue + + overlaps: list[set[int]] = [] + survivors: list[set[int]] = [] + for existing in merged: + if existing & candidate: + overlaps.append(existing) + else: + survivors.append(existing) + + for overlap in overlaps: + candidate.update(overlap) + + survivors.append(candidate) + merged = survivors + + return [sorted(group) for group in merged] + + +def build_periodic_groups( + points: np.ndarray, + bbox_min: np.ndarray, + bbox_max: np.ndarray, + periodic_axes: Iterable[bool | float], + *, + tolerance: float, +) -> list[list[int]]: + """Build periodic point groups by matching opposite boundary points.""" + + if points.size == 0: + return [] + + periodic_flags = [bool(value) for value in periodic_axes] + raw_groups: list[list[int]] = [] + + for axis, enabled in enumerate(periodic_flags): + if not enabled: + continue + + min_mask = np.isclose(points[:, axis], bbox_min[axis], atol=tolerance, rtol=0.0) + max_mask = np.isclose(points[:, axis], bbox_max[axis], atol=tolerance, rtol=0.0) + if not np.any(min_mask) or not np.any(max_mask): + continue + + other_axes = [other for other in range(points.shape[1]) if other != axis] + min_points = np.flatnonzero(min_mask) + max_points = np.flatnonzero(max_mask) + + max_lookup: defaultdict[tuple[float, ...], list[int]] = defaultdict(list) + for index in max_points: + max_lookup[_point_key(points[index, other_axes])].append(int(index)) + + for index in min_points: + key = _point_key(points[index, other_axes]) + for partner in max_lookup.get(key, []): + raw_groups.append([int(index), int(partner)]) + + return merge_periodic_groups(raw_groups) diff --git a/src/nmesh/mesher/relaxation/__init__.py b/src/nmesh/mesher/relaxation/__init__.py new file mode 100644 index 0000000..70bab4a --- /dev/null +++ b/src/nmesh/mesher/relaxation/__init__.py @@ -0,0 +1,14 @@ +"""Relaxation meshing package for the pure-Python section 4 port.""" + +from .density import _compile_density_function +from .engine import RelaxationEngine, mesh_bodies_raw +from .geometry import FemGeometry, fem_geometry_from_bodies +from .topology import assemble_raw_mesh + +__all__ = [ + "FemGeometry", + "RelaxationEngine", + "assemble_raw_mesh", + "fem_geometry_from_bodies", + "mesh_bodies_raw", +] diff --git a/src/nmesh/mesher/relaxation/_constants.py b/src/nmesh/mesher/relaxation/_constants.py new file mode 100644 index 0000000..c5b8b79 --- /dev/null +++ b/src/nmesh/mesher/relaxation/_constants.py @@ -0,0 +1,18 @@ +"""Internal constants for the relaxation meshing package.""" + +from ...utils.constants import BOUNDARY_FUZZ + +DENSITY_EPSILON = 1.0e-12 +DEFAULT_RNG_SEED = 97 +STATE_FIXED = 0 +STATE_MOBILE = 1 +STATE_SIMPLE = 2 + +__all__ = [ + "BOUNDARY_FUZZ", + "DEFAULT_RNG_SEED", + "DENSITY_EPSILON", + "STATE_FIXED", + "STATE_MOBILE", + "STATE_SIMPLE", +] diff --git a/src/nmesh/mesher/relaxation/_types.py b/src/nmesh/mesher/relaxation/_types.py new file mode 100644 index 0000000..b4a5f7f --- /dev/null +++ b/src/nmesh/mesher/relaxation/_types.py @@ -0,0 +1,14 @@ +"""Internal type aliases for the relaxation meshing package.""" + +from typing import Any, Callable, TypeAlias + +import numpy as np + +from ...utils.types import FloatArray +from ..driver import MeshEngineStatus + +DensityFunction: TypeAlias = Callable[[FloatArray], float] +RegionFunction: TypeAlias = Callable[[FloatArray], np.ndarray] +EngineResult: TypeAlias = tuple[MeshEngineStatus, Any] + +__all__ = ["DensityFunction", "EngineResult", "FloatArray", "RegionFunction"] diff --git a/src/nmesh/mesher/relaxation/density.py b/src/nmesh/mesher/relaxation/density.py new file mode 100644 index 0000000..9d98f17 --- /dev/null +++ b/src/nmesh/mesher/relaxation/density.py @@ -0,0 +1,131 @@ +"""Density-function parsing and compilation helpers.""" + +from __future__ import annotations + +import math +import re + +import numpy as np + +from ._constants import DENSITY_EPSILON +from ._types import DensityFunction, FloatArray + + +def _strip_c_comments(source: str) -> str: + """Remove C-style block comments from a legacy density script.""" + + return re.sub(r"/\*.*?\*/", "", source, flags=re.S) + + +def _translate_density_source(source: str) -> str: + """Translate a small subset of legacy C-like density syntax into Python.""" + + cleaned = _strip_c_comments(source) + cleaned = cleaned.replace("\r\n", "\n") + cleaned = cleaned.replace("{", "\n{\n").replace("}", "\n}\n").replace(";", ";\n") + lines = [line.strip() for line in cleaned.splitlines() if line.strip()] + + py_lines = ["def __compiled_density(x):", " density = 1.0"] + indent = 1 + + for line in lines: + if line == "{": + indent += 1 + continue + if line == "}": + indent = max(1, indent - 1) + continue + + translated = line.rstrip(";").strip() + translated = re.sub(r"\bdouble\s+", "", translated) + translated = translated.replace("&&", " and ").replace("||", " or ") + translated = re.sub(r"(?=!])!(?!=)", " not ", translated) + translated = re.sub(r"\bexp\s*\(", "math.exp(", translated) + translated = re.sub(r"\bsqrt\s*\(", "math.sqrt(", translated) + translated = re.sub(r"\bsin\s*\(", "math.sin(", translated) + translated = re.sub(r"\bcos\s*\(", "math.cos(", translated) + translated = re.sub(r"\btan\s*\(", "math.tan(", translated) + translated = re.sub(r"\bpow\s*\(", "math.pow(", translated) + + if translated.startswith("if "): + condition = translated[2:].strip() + if condition.startswith("(") and condition.endswith(")"): + condition = condition[1:-1].strip() + py_lines.append(f"{' ' * indent}if {condition}:") + continue + + if translated == "else": + py_lines.append(f"{' ' * indent}else:") + continue + + py_lines.append(f"{' ' * indent}{translated}") + + py_lines.append(" return float(density)") + return "\n".join(py_lines) + + +def _compile_density_function(density: str | DensityFunction | None) -> DensityFunction: + """Compile a legacy density string or callable into a normalized density function.""" + + if density is None or density == "": + return lambda _point: 1.0 + + if callable(density): + + def density_wrapper(point: FloatArray) -> float: + """Wrap a user density callable with array coercion and clamping.""" + + value = float(density(np.asarray(point, dtype=float))) + return max(value, DENSITY_EPSILON) + + return density_wrapper + + source = str(density).strip() + try: + expression = re.search(r"density\s*=\s*(.+?)\s*;", source, flags=re.S) + if expression is not None and "if" not in source and "{" not in source: + code = compile(expression.group(1), "", "eval") + + def expression_density(point: FloatArray) -> float: + """Evaluate a single-expression density snippet.""" + + x = np.asarray(point, dtype=float) + scope = { + "x": x, + "math": math, + "exp": math.exp, + "sqrt": math.sqrt, + "sin": math.sin, + "cos": math.cos, + "tan": math.tan, + "abs": abs, + "min": min, + "max": max, + } + value = float(eval(code, {"__builtins__": {}}, scope)) + return max(value, DENSITY_EPSILON) + + return expression_density + + translated = _translate_density_source(source) + namespace: dict[str, object] = { + "__builtins__": {}, + "math": math, + "float": float, + "abs": abs, + "min": min, + "max": max, + } + exec(translated, namespace, namespace) + compiled = namespace["__compiled_density"] + + def script_density(point: FloatArray) -> float: + """Evaluate a translated multi-line density script.""" + + x = np.asarray(point, dtype=float) + value = float(compiled(x)) + return max(value, DENSITY_EPSILON) + + return script_density + except Exception: + return lambda _point: 1.0 diff --git a/src/nmesh/mesher/relaxation/engine.py b/src/nmesh/mesher/relaxation/engine.py new file mode 100644 index 0000000..e218714 --- /dev/null +++ b/src/nmesh/mesher/relaxation/engine.py @@ -0,0 +1,281 @@ +"""Iterative meshing engine for the pure-Python relaxation pipeline.""" + +from __future__ import annotations + +import math +from typing import Any + +import numpy as np + +from ...backend import RawMesh +from ...geometry.primitives import Body +from ..driver import MeshEngineCommand, MeshEngineStatus +from ..meshing_parameters import PointFate, default_handle_point_density_fun +from ._constants import BOUNDARY_FUZZ, DEFAULT_RNG_SEED, STATE_MOBILE +from ._types import DensityFunction, EngineResult, FloatArray +from .geometry import FemGeometry, fem_geometry_from_bodies +from .seeding import _as_float_array, _dedupe_points, _prepare_initial_points +from .topology import _callback_mesh_info, assemble_raw_mesh + + +class RelaxationEngine: + """Small iterative meshing engine that exposes the legacy driver protocol.""" + + def __init__( + self, + geometry: FemGeometry, + mesher: dict[str, Any], + a0: float, + fixed_points: FloatArray, + mobile_points: FloatArray, + simply_points: FloatArray, + periodic: list[float] | list[bool], + *, + rng: np.random.Generator | None = None, + ) -> None: + """Initialize the engine state from geometry, seeds, and mesher parameters.""" + + self.geometry = geometry + self.params = dict(mesher.get("parameters", {})) + self.a0 = float(a0) + self.periodic = list(periodic) + self.rng = rng or np.random.default_rng(DEFAULT_RNG_SEED) + self.points, self.states = _prepare_initial_points( + geometry, + self.a0, + fixed_points, + mobile_points, + simply_points, + self.periodic, + self.rng, + ) + self.step = 0 + self.finished = False + self.last_max_displacement = math.inf + self.cached_raw_mesh = assemble_raw_mesh(self.points, self.geometry, self.periodic) + + @property + def max_steps(self) -> int: + """Return the capped step budget for this Python port.""" + + return min(int(self.params.get("controller_step_limit_max", 1000)), 60) + + @property + def tolerated_rel_move(self) -> float: + """Return the relative movement threshold used for convergence.""" + + return float(self.params.get("controller_tolerated_rel_movement", 0.002)) + + @property + def time_step_scale(self) -> float: + """Return the scale factor applied to each relaxation displacement.""" + + return float(self.params.get("controller_time_step_scale", 0.1)) + + def _rebuild_mesh(self) -> None: + """Refresh the cached ``RawMesh`` snapshot from the current point cloud.""" + + self.cached_raw_mesh = assemble_raw_mesh(self.points, self.geometry, self.periodic) + + def _neighbor_map(self) -> list[list[int]]: + """Derive point adjacency from the cached simplices.""" + + neighbors: list[set[int]] = [set() for _ in range(len(self.points))] + for simplex in self.cached_raw_mesh.simplices: + for index in simplex: + neighbors[index].update(int(item) for item in simplex if item != index) + return [sorted(group) for group in neighbors] + + def _attempt_add_delete_points(self, displacements: FloatArray) -> None: + """Apply the mesher density heuristic to add or remove mobile points.""" + + if len(self.points) == 0: + return + + handler = self.params.get("handle_point_density_fun", default_handle_point_density_fun) + thresh_add = float(self.params.get("controller_thresh_add", 1.0)) + thresh_del = float(self.params.get("controller_thresh_del", 2.0)) + neighbors = self._neighbor_map() + + additions: list[FloatArray] = [] + removals: list[int] = [] + + for index, state in enumerate(self.states): + if state != STATE_MOBILE or index in removals: + continue + + neigh = neighbors[index] + if not neigh: + continue + + point = self.points[index] + neigh_coords = self.points[np.asarray(neigh, dtype=int)] + distances = np.linalg.norm(neigh_coords - point, axis=1) + mean_distance = float(np.mean(distances)) + density_here = self.geometry.density_at(point) + ideal_distance = self.a0 / (density_here ** (1.0 / max(self.geometry.dim, 1))) + avg_density = ideal_distance / max(mean_distance, BOUNDARY_FUZZ) + avg_force = float(np.linalg.norm(displacements[index])) + + fate = handler(self.rng, (avg_density, avg_force), thresh_add, thresh_del) + if fate == PointFate.ADD_ANOTHER: + farthest = neigh_coords[int(np.argmax(distances))] + candidate = 0.5 * (point + farthest) + if self.geometry.classify_points(candidate.reshape(1, -1))[0] >= 0: + additions.append(candidate) + elif fate == PointFate.DELETE and len(self.points) - len(removals) > self.geometry.dim + 1: + removals.append(index) + + if removals: + keep = np.ones(len(self.points), dtype=bool) + keep[np.asarray(removals, dtype=int)] = False + self.points = self.points[keep] + self.states = self.states[keep] + + if additions: + additions_arr = _dedupe_points(np.asarray(additions, dtype=float)) + if len(additions_arr) > 0: + self.points = np.vstack((self.points, additions_arr)) + self.states = np.concatenate( + ( + self.states, + np.full(len(additions_arr), STATE_MOBILE, dtype=int), + ) + ) + + def _step_once(self) -> None: + """Execute one relaxation step over all currently mobile points.""" + + if len(self.points) < self.geometry.dim + 1: + self.finished = True + return + + self._rebuild_mesh() + neighbors = self._neighbor_map() + new_points = np.array(self.points, copy=True) + displacements = np.zeros_like(new_points) + max_displacement = 0.0 + + for index, state in enumerate(self.states): + if state != STATE_MOBILE: + continue + + neigh = neighbors[index] + if not neigh: + continue + + point = self.points[index] + neigh_coords = self.points[np.asarray(neigh, dtype=int)] + deltas = neigh_coords - point + distances = np.linalg.norm(deltas, axis=1) + nonzero = distances > BOUNDARY_FUZZ + if not np.any(nonzero): + continue + + deltas = deltas[nonzero] + distances = distances[nonzero] + density_here = self.geometry.density_at(point) + ideal_distance = self.a0 / (density_here ** (1.0 / max(self.geometry.dim, 1))) + scale = (distances - ideal_distance) / distances + displacement = np.mean(deltas * scale[:, None], axis=0) + displacement *= self.time_step_scale + max_norm = self.a0 * float(self.params.get("controller_movement_max_freedom", 3.0)) * 0.05 + disp_norm = float(np.linalg.norm(displacement)) + if disp_norm > max_norm > 0.0: + displacement *= max_norm / disp_norm + + candidate = point + displacement + if self.geometry.classify_points(candidate.reshape(1, -1))[0] >= 0: + new_points[index] = candidate + displacements[index] = displacement + max_displacement = max(max_displacement, float(np.linalg.norm(displacement))) + + self.points = new_points + self.last_max_displacement = max_displacement / max(self.a0, BOUNDARY_FUZZ) + + if self.step > 0 and self.step % 5 == 0: + self._attempt_add_delete_points(displacements) + + self._rebuild_mesh() + if self.last_max_displacement <= self.tolerated_rel_move: + self.finished = True + + def callback_payload(self) -> list[list[Any]]: + """Return the callback payload for the engine's current mesh snapshot.""" + + return _callback_mesh_info(self.cached_raw_mesh) + + def run(self, command: MeshEngineCommand) -> EngineResult: + """Handle one driver command and return the next engine status tuple.""" + + if command == MeshEngineCommand.DO_EXTRACT: + self._rebuild_mesh() + return MeshEngineStatus.PRODUCED_INTERMEDIATE_MESH, ( + self.callback_payload(), + self.run, + ) + + if command != MeshEngineCommand.DO_STEP: + raise ValueError(f"Unsupported mesh engine command: {command!r}") + + if self.finished or self.step >= self.max_steps: + self._rebuild_mesh() + return MeshEngineStatus.FINISHED_STEP_LIMIT_REACHED, None + + self.step += 1 + self._step_once() + + if self.finished: + self._rebuild_mesh() + return MeshEngineStatus.FINISHED_FORCE_EQUILIBRIUM_REACHED, None + + if self.step >= self.max_steps: + self._rebuild_mesh() + return MeshEngineStatus.FINISHED_STEP_LIMIT_REACHED, None + + return MeshEngineStatus.CAN_CONTINUE, self.run + + +def mesh_bodies_raw( + gendriver: Any, + mesher: dict[str, Any], + bb_min: list[float], + bb_max: list[float], + mesh_ext: int, + objects: list[Body], + a0: float, + density: str | DensityFunction | None, + fixed: list[list[float]], + mobile: list[list[float]], + simply: list[list[float]], + periodic: list[float] | list[bool], + hints: list[list[Any]], +) -> RawMesh: + """Public backend entry point for meshing bodies with the Python engine.""" + + geometry = fem_geometry_from_bodies( + (np.asarray(bb_min, dtype=float), np.asarray(bb_max, dtype=float)), + list(objects), + list(hints), + density=density, + mesh_exterior=bool(mesh_ext), + ) + + engine = RelaxationEngine( + geometry, + mesher, + float(a0), + _as_float_array(fixed, dim=geometry.dim), + _as_float_array(mobile, dim=geometry.dim), + _as_float_array(simply, dim=geometry.dim), + list(periodic), + rng=np.random.default_rng(DEFAULT_RNG_SEED), + ) + + if callable(gendriver): + gendriver(engine.run) + else: + engine.run(MeshEngineCommand.DO_STEP) + + engine._rebuild_mesh() + return engine.cached_raw_mesh diff --git a/src/nmesh/mesher/relaxation/geometry.py b/src/nmesh/mesher/relaxation/geometry.py new file mode 100644 index 0000000..d039ee3 --- /dev/null +++ b/src/nmesh/mesher/relaxation/geometry.py @@ -0,0 +1,166 @@ +"""Geometry modeling helpers for the relaxation meshing pipeline.""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any + +import numpy as np + +from ...backend import RawMesh +from ...geometry.primitives import Body +from ._constants import BOUNDARY_FUZZ, DENSITY_EPSILON +from ._types import DensityFunction, FloatArray, RegionFunction +from .density import _compile_density_function + + +@dataclass(frozen=True, slots=True) +class FemGeometry: + """Geometry bundle used by the Python meshing engine.""" + + dim: int + bbox_min: FloatArray + bbox_max: FloatArray + density_fun: DensityFunction + region_ids: tuple[int, ...] + region_functions: tuple[RegionFunction, ...] + piece_hints: tuple[FloatArray, ...] + mesh_exterior: bool + + def points_in_bbox(self, points: FloatArray) -> np.ndarray: + """Return a mask showing which points lie inside the bounding box.""" + + return np.all( + (points >= self.bbox_min - BOUNDARY_FUZZ) + & (points <= self.bbox_max + BOUNDARY_FUZZ), + axis=1, + ) + + def classify_points(self, points: FloatArray) -> np.ndarray: + """Return the first matching region id for each point, or -1 outside.""" + + if points.size == 0: + return np.empty(0, dtype=int) + + coords = np.asarray(points, dtype=float) + assigned = np.full(len(coords), -1, dtype=int) + bbox_mask = self.points_in_bbox(coords) + + for region_id, region_fun in zip(self.region_ids, self.region_functions): + mask = bbox_mask & region_fun(coords) & (assigned < 0) + assigned[mask] = region_id + + return assigned + + def density_at(self, point: FloatArray) -> float: + """Evaluate the density function with a small positive lower bound.""" + + return max(float(self.density_fun(point)), DENSITY_EPSILON) + + +def _raw_mesh_points(raw_mesh: RawMesh, dim: int) -> FloatArray: + """Return mesh points as a validated floating-point array.""" + + coords = np.asarray(raw_mesh.points, dtype=float) + if coords.size == 0: + return np.empty((0, dim), dtype=float) + if coords.ndim == 1: + coords = coords.reshape(1, dim) + if coords.shape[1] != dim: + raise ValueError(f"Expected hint mesh points with dimension {dim}, got {coords.shape[1]}") + return coords.astype(float, copy=False) + + +def _dedupe_hint_points(points: FloatArray) -> FloatArray: + """Drop duplicate hint points while preserving first-seen order.""" + + if len(points) == 0: + return points + + keep_indices: list[int] = [] + seen: set[tuple[float, ...]] = set() + for index, point in enumerate(points): + key = tuple(np.round(np.asarray(point, dtype=float), decimals=10).tolist()) + if key in seen: + continue + seen.add(key) + keep_indices.append(index) + return points[np.asarray(keep_indices, dtype=int)] + + +def fem_geometry_from_bodies( + bounding_box: tuple[FloatArray, FloatArray], + bodies: list[Body], + hints: list[list[Any]], + *, + density: str | DensityFunction | None = None, + mesh_exterior: bool = False, +) -> FemGeometry: + """Construct the geometry bundle consumed by the Python meshing engine.""" + + bbox_min = np.asarray(bounding_box[0], dtype=float) + bbox_max = np.asarray(bounding_box[1], dtype=float) + dim = int(len(bbox_min)) + density_fun = _compile_density_function(density) + + body_functions: list[RegionFunction] = [] + piece_hints: list[FloatArray] = [] + region_ids: list[int] = [] + next_region_id = 1 + + for body in bodies: + body_functions.append( + lambda points, member=body: np.asarray(member.evaluate(points), dtype=float) + >= -BOUNDARY_FUZZ + ) + piece_hints.append(np.empty((0, dim), dtype=float)) + region_ids.append(next_region_id) + next_region_id += 1 + + for hint_mesh, hint_body in hints: + hint_points = _raw_mesh_points(hint_mesh, dim) + if len(hint_points) > 0: + mask = np.asarray(hint_body.evaluate(hint_points), dtype=float) >= -BOUNDARY_FUZZ + hint_points = _dedupe_hint_points(hint_points[mask]) + body_functions.append( + lambda points, member=hint_body: np.asarray(member.evaluate(points), dtype=float) + >= -BOUNDARY_FUZZ + ) + piece_hints.append(hint_points) + region_ids.append(next_region_id) + next_region_id += 1 + + region_functions: list[RegionFunction] = [] + ordered_region_ids: list[int] = [] + + if mesh_exterior: + + def exterior(points: FloatArray) -> np.ndarray: + """Classify points that lie inside the box but outside all bodies.""" + + bbox_mask = np.all( + (points >= bbox_min - BOUNDARY_FUZZ) + & (points <= bbox_max + BOUNDARY_FUZZ), + axis=1, + ) + inside_any = np.zeros(len(points), dtype=bool) + for body_fun in body_functions: + inside_any |= body_fun(points) + return bbox_mask & ~inside_any + + region_functions.append(exterior) + ordered_region_ids.append(0) + + region_functions.extend(body_functions) + ordered_region_ids.extend(region_ids) + + return FemGeometry( + dim=dim, + bbox_min=bbox_min, + bbox_max=bbox_max, + density_fun=density_fun, + region_ids=tuple(ordered_region_ids), + region_functions=tuple(region_functions), + piece_hints=tuple(piece_hints), + mesh_exterior=bool(mesh_exterior), + ) diff --git a/src/nmesh/mesher/relaxation/seeding.py b/src/nmesh/mesher/relaxation/seeding.py new file mode 100644 index 0000000..17ee460 --- /dev/null +++ b/src/nmesh/mesher/relaxation/seeding.py @@ -0,0 +1,225 @@ +"""Seed-point preparation helpers for the relaxation meshing pipeline.""" + +from __future__ import annotations + +from typing import Any + +import numpy as np + +from ._constants import BOUNDARY_FUZZ, DEFAULT_RNG_SEED, DENSITY_EPSILON, STATE_FIXED, STATE_MOBILE, STATE_SIMPLE +from ._types import FloatArray +from .geometry import FemGeometry + + +def _as_float_array(points: Any, dim: int | None = None) -> FloatArray: + """Convert point-like input into a 2D float array with optional dimension validation.""" + + if points is None: + target_dim = 0 if dim is None else dim + return np.empty((0, target_dim), dtype=float) + + coords = np.asarray(points, dtype=float) + if coords.size == 0: + target_dim = 0 if dim is None else dim + return np.empty((0, target_dim), dtype=float) + + if coords.ndim == 1: + if dim is None: + dim = int(coords.shape[0]) + coords = coords.reshape(1, dim) + + if dim is not None and coords.shape[1] != dim: + raise ValueError(f"Expected points with dimension {dim}, got {coords.shape[1]}") + + return coords.astype(float, copy=False) + + +def _point_key(point: FloatArray, decimals: int = 10) -> tuple[float, ...]: + """Return a rounded tuple key for deduplication of point coordinates.""" + + return tuple(np.round(np.asarray(point, dtype=float), decimals=decimals).tolist()) + + +def _dedupe_points(points: FloatArray) -> FloatArray: + """Drop duplicate points while preserving first-seen order.""" + + if len(points) == 0: + return points + + keep_indices: list[int] = [] + seen: set[tuple[float, ...]] = set() + for index, point in enumerate(points): + key = _point_key(point) + if key in seen: + continue + seen.add(key) + keep_indices.append(index) + return points[np.asarray(keep_indices, dtype=int)] + + +def _dedupe_fixed_mobile( + fixed_points: FloatArray, + mobile_points: FloatArray, + simply_points: FloatArray, +) -> tuple[FloatArray, FloatArray, FloatArray]: + """Deduplicate seed points across fixed, mobile, and simply categories.""" + + seen: set[tuple[float, ...]] = set() + + def filter_points(points: FloatArray) -> FloatArray: + """Return only points not seen in an earlier category.""" + + keep_indices: list[int] = [] + for index, point in enumerate(points): + key = _point_key(point) + if key in seen: + continue + seen.add(key) + keep_indices.append(index) + if not keep_indices: + return np.empty((0, points.shape[1]), dtype=float) + return points[np.asarray(keep_indices, dtype=int)] + + return ( + filter_points(fixed_points), + filter_points(mobile_points), + filter_points(simply_points), + ) + + +def _grid_candidate_points(geometry: FemGeometry, a0: float) -> FloatArray: + """Generate deterministic candidate seed points over the bounding box.""" + + extents = np.maximum(geometry.bbox_max - geometry.bbox_min, a0) + sample_points = [geometry.bbox_min, geometry.bbox_max, 0.5 * (geometry.bbox_min + geometry.bbox_max)] + max_density = max(geometry.density_at(point) for point in sample_points) + density_scale = min(max_density ** (1.0 / max(geometry.dim, 1)), 2.0) + step = max(a0 / density_scale, a0 * 0.5) + + counts = np.maximum(np.floor(extents / step).astype(int) + 1, 2) + candidate_count = int(np.prod(counts)) + if candidate_count > 15_000: + approximate = min(8_000, max(500, int(np.prod(extents / max(a0, DENSITY_EPSILON)) * 4))) + rng = np.random.default_rng(DEFAULT_RNG_SEED) + return rng.uniform(geometry.bbox_min, geometry.bbox_max, size=(approximate, geometry.dim)) + + axes = [ + np.linspace(geometry.bbox_min[axis], geometry.bbox_max[axis], int(counts[axis])) + for axis in range(geometry.dim) + ] + grid = np.meshgrid(*axes, indexing="ij") + return np.stack(grid, axis=-1).reshape(-1, geometry.dim) + + +def _keep_point_for_density(point: FloatArray, density_here: float, rng: np.random.Generator) -> bool: + """Return whether a candidate point survives density-based thinning.""" + + _ = point + if density_here >= 1.0: + return True + return bool(rng.random() <= density_here) + + +def _prepare_initial_points( + geometry: FemGeometry, + a0: float, + fixed_points: FloatArray, + mobile_points: FloatArray, + simply_points: FloatArray, + periodic: list[float] | list[bool], + rng: np.random.Generator, +) -> tuple[FloatArray, np.ndarray]: + """Prepare the initial point cloud and point-state array for relaxation.""" + + fixed_points, mobile_points, simply_points = _dedupe_fixed_mobile( + fixed_points, mobile_points, simply_points + ) + + def filter_relevant(points: FloatArray) -> FloatArray: + """Keep only points that belong to a meshed region.""" + + if len(points) == 0: + return points + mask = geometry.classify_points(points) >= 0 + return points[mask] + + fixed_points = filter_relevant(fixed_points) + mobile_points = filter_relevant(mobile_points) + simply_points = filter_relevant(simply_points) + + candidates = _grid_candidate_points(geometry, a0) + region_ids = geometry.classify_points(candidates) + candidates = candidates[region_ids >= 0] + + if len(candidates) > 0: + candidate_keys = {_point_key(point) for point in fixed_points} + candidate_keys.update(_point_key(point) for point in mobile_points) + candidate_keys.update(_point_key(point) for point in simply_points) + + selected_candidates: list[FloatArray] = [] + for point in candidates: + key = _point_key(point) + if key in candidate_keys: + continue + density_here = geometry.density_at(point) + if not _keep_point_for_density(point, density_here, rng): + continue + selected_candidates.append(point) + candidate_keys.add(key) + generated_points = ( + np.asarray(selected_candidates, dtype=float) + if selected_candidates + else np.empty((0, geometry.dim), dtype=float) + ) + else: + generated_points = np.empty((0, geometry.dim), dtype=float) + + hint_points = [ + hint_points + for hint_points in geometry.piece_hints + if len(hint_points) > 0 + ] + if hint_points: + hint_block = _dedupe_points(np.vstack(hint_points)) + hint_block = filter_relevant(hint_block) + else: + hint_block = np.empty((0, geometry.dim), dtype=float) + + all_points = np.vstack( + ( + fixed_points, + simply_points, + mobile_points, + hint_block, + generated_points, + ) + ) + + states = np.concatenate( + ( + np.full(len(fixed_points), STATE_FIXED, dtype=int), + np.full(len(simply_points), STATE_SIMPLE, dtype=int), + np.full(len(mobile_points), STATE_MOBILE, dtype=int), + np.full(len(hint_block), STATE_FIXED, dtype=int), + np.full(len(generated_points), STATE_MOBILE, dtype=int), + ) + ) + + if len(all_points) == 0: + return all_points, states + + periodic_flags = [bool(value) for value in periodic] + periodic_mask = np.zeros(len(all_points), dtype=bool) + tolerance = max(0.05 * a0, BOUNDARY_FUZZ) + for axis, enabled in enumerate(periodic_flags): + if not enabled: + continue + periodic_mask |= np.isclose( + all_points[:, axis], geometry.bbox_min[axis], atol=tolerance, rtol=0.0 + ) + periodic_mask |= np.isclose( + all_points[:, axis], geometry.bbox_max[axis], atol=tolerance, rtol=0.0 + ) + states[periodic_mask] = STATE_FIXED + + return all_points, states diff --git a/src/nmesh/mesher/relaxation/topology.py b/src/nmesh/mesher/relaxation/topology.py new file mode 100644 index 0000000..7772a02 --- /dev/null +++ b/src/nmesh/mesher/relaxation/topology.py @@ -0,0 +1,193 @@ +"""Topology extraction and RawMesh assembly helpers.""" + +from __future__ import annotations + +import math +from itertools import combinations +from typing import Any + +import numpy as np +from scipy.spatial import Delaunay, QhullError + +from ...backend import RawMesh +from ..periodic import build_periodic_groups +from ._constants import BOUNDARY_FUZZ, STATE_MOBILE +from ._types import FloatArray +from .geometry import FemGeometry + + +def _simplex_measures(points: FloatArray, simplices: np.ndarray, dim: int) -> FloatArray: + """Compute 1D lengths or higher-dimensional simplex volumes.""" + + if len(simplices) == 0: + return np.empty(0, dtype=float) + + if dim == 1: + return np.abs(points[simplices[:, 1], 0] - points[simplices[:, 0], 0]) + + matrices = points[simplices[:, 1:]] - points[simplices[:, [0]]] + determinants = np.linalg.det(matrices) + return np.abs(determinants) / math.factorial(dim) + + +def _surface_faces(simplices: np.ndarray, regions: np.ndarray, dim: int) -> list[list[int]]: + """Extract boundary faces by counting face ownership across simplices.""" + + if len(simplices) == 0: + return [] + + face_map: dict[tuple[int, ...], set[int]] = {} + face_size = dim + for simplex, region in zip(simplices, regions): + for face in combinations(simplex.tolist(), face_size): + key = tuple(sorted(face)) + face_map.setdefault(key, set()).add(int(region)) + + counts: dict[tuple[int, ...], int] = {} + for simplex in simplices: + for face in combinations(simplex.tolist(), face_size): + key = tuple(sorted(face)) + counts[key] = counts.get(key, 0) + 1 + + surfaces: list[list[int]] = [] + for face, count in counts.items(): + if count == 1 or len(face_map.get(face, set())) > 1: + surfaces.append(list(face)) + return surfaces + + +def _unique_links(simplices: np.ndarray) -> list[tuple[int, int]]: + """Return the sorted unique undirected edges induced by the simplices.""" + + links: set[tuple[int, int]] = set() + for simplex in simplices: + for edge_start, edge_end in combinations(simplex.tolist(), 2): + start = int(edge_start) + end = int(edge_end) + if start <= end: + links.add((start, end)) + else: + links.add((end, start)) + return sorted(links) + + +def _point_regions(point_count: int, simplices: np.ndarray, regions: np.ndarray) -> list[list[int]]: + """Build the region-membership list for each mesh point.""" + + memberships: list[set[int]] = [set() for _ in range(point_count)] + for simplex, region in zip(simplices, regions): + for point_index in simplex: + memberships[int(point_index)].add(int(region)) + return [sorted(group) for group in memberships] + + +def _region_volumes(region_ids: np.ndarray, measures: FloatArray) -> list[float]: + """Aggregate simplex measures into per-region total volumes.""" + + if len(region_ids) == 0: + return [] + + order = sorted({int(region) for region in region_ids}) + totals = {region: 0.0 for region in order} + for region, measure in zip(region_ids, measures): + totals[int(region)] += float(measure) + return [totals[region] for region in order] + + +def _triangulate_points(points: FloatArray, dim: int, states: np.ndarray | None = None) -> np.ndarray: + """Triangulate the point cloud, retrying with light jitter for degenerate inputs.""" + + if len(points) < dim + 1: + return np.empty((0, dim + 1), dtype=int) + + if dim == 1: + order = np.argsort(points[:, 0], kind="mergesort") + return np.column_stack((order[:-1], order[1:])).astype(int) + + try: + return Delaunay(points).simplices.astype(int, copy=False) + except QhullError: + jittered = np.array(points, copy=True) + movable = np.ones(len(points), dtype=bool) if states is None else states == STATE_MOBILE + if np.any(movable): + amplitudes = np.linspace(1.0e-9, 1.0e-8, np.count_nonzero(movable)) + offsets = np.column_stack( + [ + amplitudes * np.sin(np.arange(1, len(amplitudes) + 1) * (axis + 1)) + for axis in range(dim) + ] + ) + jittered[movable] += offsets + try: + return Delaunay(jittered).simplices.astype(int, copy=False) + except QhullError: + return np.empty((0, dim + 1), dtype=int) + + +def assemble_raw_mesh( + points: FloatArray, + geometry: FemGeometry, + periodic: list[float] | list[bool], +) -> RawMesh: + """Assemble a ``RawMesh`` from relaxed points and geometry classification.""" + + coords = np.asarray(points, dtype=float) + dim = geometry.dim + simplices = _triangulate_points(coords, dim) + + if len(simplices) > 0: + centroids = np.mean(coords[simplices], axis=1) + region_ids = geometry.classify_points(centroids) + measures = _simplex_measures(coords, simplices, dim) + keep = (region_ids >= 0) & (measures > BOUNDARY_FUZZ) + simplices = simplices[keep] + region_ids = region_ids[keep] + measures = measures[keep] + else: + region_ids = np.empty(0, dtype=int) + measures = np.empty(0, dtype=float) + + surfaces = _surface_faces(simplices, region_ids, dim) + links = _unique_links(simplices) + point_regions = _point_regions(len(coords), simplices, region_ids) + region_volumes = _region_volumes(region_ids, measures) + periodic_groups = build_periodic_groups( + coords, + geometry.bbox_min, + geometry.bbox_max, + periodic, + tolerance=BOUNDARY_FUZZ, + ) + + return RawMesh( + points=coords.tolist(), + simplices=simplices.tolist(), + regions=region_ids.astype(int).tolist(), + point_regions=point_regions, + surfaces=surfaces, + links=links, + region_volumes=region_volumes, + periodic_point_indices=periodic_groups, + permutation=list(range(len(coords))), + dim=dim, + ) + + +def _callback_mesh_info(raw_mesh: RawMesh) -> list[list[Any]]: + """Build the legacy-style callback payload expected by old consumers.""" + + simplices = [ + [simplex, (([], 0.0), ([], 0.0), region)] + for simplex, region in zip(raw_mesh.simplices, raw_mesh.regions) + ] + surfaces = [ + [surface, (([], 0.0), ([], 0.0), 1)] + for surface in raw_mesh.surfaces + ] + return [ + ["COORDS", "Node coordinates", raw_mesh.points], + ["LINKS", "Mesh links", raw_mesh.links], + ["POINT-BODIES", "Point region memberships", raw_mesh.point_regions], + ["SIMPLICES", "Simplex connectivity", simplices], + ["SURFACES", "Surface connectivity", surfaces], + ] diff --git a/src/nmesh/utils/constants.py b/src/nmesh/utils/constants.py index 7d34416..d6800fc 100644 --- a/src/nmesh/utils/constants.py +++ b/src/nmesh/utils/constants.py @@ -7,10 +7,14 @@ """ __all__ = [ + "BOUNDARY_FUZZ", "MIN_ABS_SCALE_FACTOR", "MIN_DIVISION_MAGNITUDE", ] +# Default tolerance used when testing whether points lie on or near boundaries. +BOUNDARY_FUZZ = 1e-6 + # Smallest allowed absolute scale factor for affine transformations. MIN_ABS_SCALE_FACTOR = 1e-12 diff --git a/tests/nmesh/test_meshing_engine.py b/tests/nmesh/test_meshing_engine.py new file mode 100644 index 0000000..48ad9ee --- /dev/null +++ b/tests/nmesh/test_meshing_engine.py @@ -0,0 +1,46 @@ +import unittest + +import numpy as np + +from nmesh import Box, Mesh +from nmesh.mesher.relaxation import _compile_density_function + + +class TestMeshingEngine(unittest.TestCase): + def test_density_string_translation_supports_c_style_blocks(self): + density = _compile_density_function( + """ + double upper=7.1; + double lower=5.9; + if ((x[0] < 2.1) && (x[0] > -2.1) && (x[1] > lower) && (x[1] < upper)) + {density = 4.0;} + else { + double sigma=1.0; + double xdev = 0.0-x[0]; + double ydev = 0.0-x[1]; + double rdev2 = xdev*xdev+ydev*ydev; + density=1.0+4.0*exp(-rdev2/(sigma*sigma)); + } + """ + ) + + self.assertAlmostEqual(density(np.asarray([0.0, 6.0])), 4.0) + self.assertAlmostEqual(density(np.asarray([0.0, 0.0])), 5.0) + + def test_meshing_is_deterministic_for_identical_inputs(self): + kwargs = dict( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[Box([0.1, 0.1], [0.9, 0.9])], + a0=0.5, + ) + + mesh_a = Mesh(**kwargs) + mesh_b = Mesh(**kwargs) + + self.assertEqual(mesh_a.points, mesh_b.points) + self.assertEqual(mesh_a.simplices, mesh_b.simplices) + self.assertEqual(mesh_a.regions, mesh_b.regions) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/nmesh_test.py b/tests/nmesh_test.py index ddf68c1..5b0858f 100644 --- a/tests/nmesh_test.py +++ b/tests/nmesh_test.py @@ -47,18 +47,72 @@ def test_csg_operations(self): d = nmesh.difference(b1, [b2]) self.assertEqual(d.dim, 3) - def test_mesh_generation_stub(self): - """Test Mesh class initialization with stubs.""" + def test_mesh_generation_produces_topology(self): + """Test Mesh class initialization now produces a populated mesh.""" bb = [[0,0,0], [1,1,1]] obj = nmesh.Box([0.2,0.2,0.2], [0.8,0.8,0.8]) m = nmesh.Mesh(bounding_box=bb, objects=[obj], a0=0.1) - self.assertEqual(str(m), "Mesh with 0 points and 0 simplices") # From stubs - - # Test properties (should return empty lists from stubs) - self.assertEqual(m.points, []) - self.assertEqual(m.simplices, []) - self.assertEqual(m.regions, []) + self.assertGreater(len(m.points), 0) + self.assertGreater(len(m.simplices), 0) + self.assertTrue(all(region == 1 for region in m.regions)) + + def test_periodic_mesh_generation(self): + """Periodic meshes should record point equivalence groups.""" + bb = [[0.0, 0.0], [1.0, 1.0]] + obj = nmesh.Box([0.25, 0.25], [0.75, 0.75]) + + m = nmesh.Mesh( + bounding_box=bb, + objects=[obj], + a0=0.5, + periodic=[True, False], + mesh_bounding_box=True, + ) + + self.assertGreater(len(m.points), 0) + self.assertGreater(len(m.simplices), 0) + self.assertGreater(len(m.periodic_point_indices), 0) + self.assertTrue(all(len(group) >= 2 for group in m.periodic_point_indices)) + + def test_hint_mesh_seeds_generation(self): + """Hints should contribute seed points to the generated mesh.""" + seed = nmesh.mesh_from_points_and_simplices( + points=[[0.2, 0.2], [0.8, 0.2], [0.5, 0.7]], + simplices_indices=[[0, 1, 2]], + simplices_regions=[3], + ) + + mesh = nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[], + mesh_bounding_box=True, + hints=[(seed, nmesh.Box([0.1, 0.1], [0.9, 0.9]))], + a0=0.5, + ) + + self.assertIn([0.2, 0.2], mesh.points) + self.assertIn([0.8, 0.2], mesh.points) + self.assertIn([0.5, 0.7], mesh.points) + + def test_callback_receives_mesh_info_payload(self): + """Callbacks should receive the legacy-style mesh_info bundle.""" + calls = [] + + def callback(piece, step, mesh_info): + calls.append((piece, step, [entry[0] for entry in mesh_info])) + + mesh = nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[nmesh.Box([0.0, 0.0], [1.0, 1.0])], + a0=0.5, + callback=(callback, 2), + ) + + self.assertGreater(len(mesh.points), 0) + self.assertGreaterEqual(len(calls), 1) + self.assertEqual(calls[0][0], 0) + self.assertEqual(calls[0][2], ["COORDS", "LINKS", "POINT-BODIES", "SIMPLICES", "SURFACES"]) def test_1d_mesh_generation(self): """Test 1D mesh generation logic.""" From 914968fd6cb2256fccad656d9d024a6dafe3b191 Mon Sep 17 00:00:00 2001 From: Emanuel Pituch <32016786+epituch@users.noreply.github.com> Date: Wed, 22 Apr 2026 22:04:52 -0700 Subject: [PATCH 2/3] Section 4 complete --- pyproject.toml | 2 +- src/nmesh/mesher/relaxation/_constants.py | 4 +- src/nmesh/mesher/relaxation/density.py | 139 ++-- src/nmesh/mesher/relaxation/engine.py | 253 +++++-- src/nmesh/mesher/relaxation/forces.py | 777 ++++++++++++++++++++++ src/nmesh/mesher/relaxation/geometry.py | 179 ++++- src/nmesh/mesher/relaxation/jit.py | 102 +++ src/nmesh/mesher/relaxation/seeding.py | 185 ++++-- src/nmesh/mesher/relaxation/topology.py | 74 ++- tests/nmesh/test_integration.py | 58 ++ tests/nmesh/test_meshing_engine.py | 165 ++++- 11 files changed, 1732 insertions(+), 206 deletions(-) create mode 100644 src/nmesh/mesher/relaxation/forces.py create mode 100644 src/nmesh/mesher/relaxation/jit.py create mode 100644 tests/nmesh/test_integration.py diff --git a/pyproject.toml b/pyproject.toml index 348a1b8..428f70b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,7 +12,7 @@ requires-python = ">=3.10" classifiers = [ "Programming Language :: Python :: 3", "Operating System :: OS Independent" ] -dependencies = ["pint", "numpy", "scipy", "tabulate", "meshio", "h5py"] +dependencies = ["pint", "numpy", "scipy", "numba", "tabulate", "meshio", "h5py"] [project.optional-dependencies] test = ["pytest", "pytest-cov", "pytest-watch"] diff --git a/src/nmesh/mesher/relaxation/_constants.py b/src/nmesh/mesher/relaxation/_constants.py index c5b8b79..3266e14 100644 --- a/src/nmesh/mesher/relaxation/_constants.py +++ b/src/nmesh/mesher/relaxation/_constants.py @@ -6,12 +6,14 @@ DEFAULT_RNG_SEED = 97 STATE_FIXED = 0 STATE_MOBILE = 1 -STATE_SIMPLE = 2 +STATE_BOUNDARY = 2 +STATE_SIMPLE = 3 __all__ = [ "BOUNDARY_FUZZ", "DEFAULT_RNG_SEED", "DENSITY_EPSILON", + "STATE_BOUNDARY", "STATE_FIXED", "STATE_MOBILE", "STATE_SIMPLE", diff --git a/src/nmesh/mesher/relaxation/density.py b/src/nmesh/mesher/relaxation/density.py index 9d98f17..6f9d4bf 100644 --- a/src/nmesh/mesher/relaxation/density.py +++ b/src/nmesh/mesher/relaxation/density.py @@ -64,6 +64,89 @@ def _translate_density_source(source: str) -> str: return "\n".join(py_lines) +def _clamp_density_value(value: float) -> float: + """Clamp a density sample to the smallest supported positive value.""" + + return max(float(value), DENSITY_EPSILON) + + +def _expression_density_scope(point: FloatArray) -> dict[str, object]: + """Build the safe eval scope for single-expression density snippets.""" + + x = np.asarray(point, dtype=float) + return { + "x": x, + "math": math, + "exp": math.exp, + "sqrt": math.sqrt, + "sin": math.sin, + "cos": math.cos, + "tan": math.tan, + "abs": abs, + "min": min, + "max": max, + } + + +def _build_callable_density(density: DensityFunction) -> DensityFunction: + """Wrap a user callable with array coercion and positive clamping.""" + + def density_wrapper(point: FloatArray) -> float: + """Evaluate a user density callable safely.""" + + return _clamp_density_value(float(density(np.asarray(point, dtype=float)))) + + return density_wrapper + + +def _compile_expression_density(source: str) -> DensityFunction | None: + """Compile a single-expression density assignment when available.""" + + expression = re.search(r"density\s*=\s*(.+?)\s*;", source, flags=re.S) + if expression is None or "if" in source or "{" in source: + return None + + code = compile(expression.group(1), "", "eval") + + def expression_density(point: FloatArray) -> float: + """Evaluate a single-expression density snippet.""" + + value = float(eval(code, {"__builtins__": {}}, _expression_density_scope(point))) + return _clamp_density_value(value) + + return expression_density + + +def _script_density_namespace() -> dict[str, object]: + """Return the restricted execution namespace for translated density scripts.""" + + return { + "__builtins__": {}, + "math": math, + "float": float, + "abs": abs, + "min": min, + "max": max, + } + + +def _compile_script_density(source: str) -> DensityFunction: + """Compile a translated multi-line legacy density script.""" + + translated = _translate_density_source(source) + namespace = _script_density_namespace() + exec(translated, namespace, namespace) + compiled = namespace["__compiled_density"] + + def script_density(point: FloatArray) -> float: + """Evaluate a translated multi-line density script.""" + + x = np.asarray(point, dtype=float) + return _clamp_density_value(float(compiled(x))) + + return script_density + + def _compile_density_function(density: str | DensityFunction | None) -> DensityFunction: """Compile a legacy density string or callable into a normalized density function.""" @@ -71,61 +154,13 @@ def _compile_density_function(density: str | DensityFunction | None) -> DensityF return lambda _point: 1.0 if callable(density): - - def density_wrapper(point: FloatArray) -> float: - """Wrap a user density callable with array coercion and clamping.""" - - value = float(density(np.asarray(point, dtype=float))) - return max(value, DENSITY_EPSILON) - - return density_wrapper + return _build_callable_density(density) source = str(density).strip() try: - expression = re.search(r"density\s*=\s*(.+?)\s*;", source, flags=re.S) - if expression is not None and "if" not in source and "{" not in source: - code = compile(expression.group(1), "", "eval") - - def expression_density(point: FloatArray) -> float: - """Evaluate a single-expression density snippet.""" - - x = np.asarray(point, dtype=float) - scope = { - "x": x, - "math": math, - "exp": math.exp, - "sqrt": math.sqrt, - "sin": math.sin, - "cos": math.cos, - "tan": math.tan, - "abs": abs, - "min": min, - "max": max, - } - value = float(eval(code, {"__builtins__": {}}, scope)) - return max(value, DENSITY_EPSILON) - + expression_density = _compile_expression_density(source) + if expression_density is not None: return expression_density - - translated = _translate_density_source(source) - namespace: dict[str, object] = { - "__builtins__": {}, - "math": math, - "float": float, - "abs": abs, - "min": min, - "max": max, - } - exec(translated, namespace, namespace) - compiled = namespace["__compiled_density"] - - def script_density(point: FloatArray) -> float: - """Evaluate a translated multi-line density script.""" - - x = np.asarray(point, dtype=float) - value = float(compiled(x)) - return max(value, DENSITY_EPSILON) - - return script_density + return _compile_script_density(source) except Exception: return lambda _point: 1.0 diff --git a/src/nmesh/mesher/relaxation/engine.py b/src/nmesh/mesher/relaxation/engine.py index e218714..e468f58 100644 --- a/src/nmesh/mesher/relaxation/engine.py +++ b/src/nmesh/mesher/relaxation/engine.py @@ -2,6 +2,7 @@ from __future__ import annotations +import logging import math from typing import Any @@ -11,12 +12,15 @@ from ...geometry.primitives import Body from ..driver import MeshEngineCommand, MeshEngineStatus from ..meshing_parameters import PointFate, default_handle_point_density_fun -from ._constants import BOUNDARY_FUZZ, DEFAULT_RNG_SEED, STATE_MOBILE +from ._constants import BOUNDARY_FUZZ, DEFAULT_RNG_SEED, STATE_BOUNDARY, STATE_FIXED, STATE_MOBILE, STATE_SIMPLE from ._types import DensityFunction, EngineResult, FloatArray +from .forces import ForceSummary, compute_forces from .geometry import FemGeometry, fem_geometry_from_bodies -from .seeding import _as_float_array, _dedupe_points, _prepare_initial_points +from .seeding import _as_float_array, _classify_dynamic_states, _dedupe_points, _prepare_initial_points from .topology import _callback_mesh_info, assemble_raw_mesh +log = logging.getLogger(__name__) + class RelaxationEngine: """Small iterative meshing engine that exposes the legacy driver protocol.""" @@ -52,7 +56,15 @@ def __init__( self.step = 0 self.finished = False self.last_max_displacement = math.inf - self.cached_raw_mesh = assemble_raw_mesh(self.points, self.geometry, self.periodic) + self.last_max_effective_force = math.inf + self._refresh_boundary_states() + self.cached_raw_mesh = assemble_raw_mesh( + self.points, + self.geometry, + self.periodic, + states=self.states, + params=self.params, + ) @property def max_steps(self) -> int: @@ -72,50 +84,103 @@ def time_step_scale(self) -> float: return float(self.params.get("controller_time_step_scale", 0.1)) + @property + def max_time_step(self) -> float: + """Return the largest time-step magnitude allowed by the controller.""" + + return float(self.params.get("controller_max_time_step", 10.0)) + + @property + def boundary_drift_tolerance(self) -> float: + """Return the tolerated increase in boundary distance for boundary points. + + We allow a small positive tolerance here rather than requiring a strict + monotonic decrease on every step. That matches the practical numerical + behavior of the Python port better and avoids jitter from floating-point + noise when a point is already very close to the surface. + """ + + return max(BOUNDARY_FUZZ, 0.02 * self.a0) + + @property + def min_equilibrium_steps(self) -> int: + """Return the minimum number of steps before equilibrium can stop the loop.""" + + configured = int(self.params.get("controller_step_limit_min", 500)) + return min(configured, min(4, self.max_steps)) + def _rebuild_mesh(self) -> None: """Refresh the cached ``RawMesh`` snapshot from the current point cloud.""" - self.cached_raw_mesh = assemble_raw_mesh(self.points, self.geometry, self.periodic) + self.cached_raw_mesh = assemble_raw_mesh( + self.points, + self.geometry, + self.periodic, + states=self.states, + params=self.params, + ) - def _neighbor_map(self) -> list[list[int]]: - """Derive point adjacency from the cached simplices.""" + def _refresh_boundary_states(self) -> None: + """Refresh dynamic points between mobile and boundary states.""" - neighbors: list[set[int]] = [set() for _ in range(len(self.points))] - for simplex in self.cached_raw_mesh.simplices: - for index in simplex: - neighbors[index].update(int(item) for item in simplex if item != index) - return [sorted(group) for group in neighbors] + dynamic_mask = np.isin(self.states, [STATE_MOBILE, STATE_BOUNDARY]) + if not np.any(dynamic_mask): + return + refreshed = _classify_dynamic_states(self.geometry, self.points[dynamic_mask], self.a0) + self.states[dynamic_mask] = refreshed + + def _effective_time_step(self, force_summary: ForceSummary) -> float: + """Compute the controller-like time step for the current force field.""" + + settling_steps = int(self.params.get("controller_initial_settling_steps", 100)) + weight_fun = self.params.get("initial_relaxation_weight_fun") + if callable(weight_fun): + relaxation_weight = float(weight_fun(self.step, settling_steps, 0.0, 1.0)) + else: + relaxation_weight = min(1.0, float(self.step) / max(float(settling_steps), 1.0)) + capped_max_time_step = relaxation_weight * self.max_time_step + if force_summary.max_effective_force <= BOUNDARY_FUZZ: + return capped_max_time_step + return min( + capped_max_time_step, + (relaxation_weight * self.time_step_scale) / force_summary.max_effective_force, + ) - def _attempt_add_delete_points(self, displacements: FloatArray) -> None: + def _attempt_add_delete_points(self, force_summary: ForceSummary) -> None: """Apply the mesher density heuristic to add or remove mobile points.""" if len(self.points) == 0: return + additions, removals = self._evaluate_point_densities(force_summary) + self._apply_point_changes(additions, removals) + self._refresh_boundary_states() + + def _evaluate_point_densities( + self, + force_summary: ForceSummary, + ) -> tuple[list[FloatArray], list[int]]: + """Decide which points to add or remove from the current cloud.""" + handler = self.params.get("handle_point_density_fun", default_handle_point_density_fun) thresh_add = float(self.params.get("controller_thresh_add", 1.0)) thresh_del = float(self.params.get("controller_thresh_del", 2.0)) - neighbors = self._neighbor_map() - additions: list[FloatArray] = [] removals: list[int] = [] for index, state in enumerate(self.states): - if state != STATE_MOBILE or index in removals: + if state not in (STATE_MOBILE, STATE_BOUNDARY) or index in removals: continue - neigh = neighbors[index] + neigh = force_summary.neighbor_map[index] if not neigh: continue point = self.points[index] neigh_coords = self.points[np.asarray(neigh, dtype=int)] distances = np.linalg.norm(neigh_coords - point, axis=1) - mean_distance = float(np.mean(distances)) - density_here = self.geometry.density_at(point) - ideal_distance = self.a0 / (density_here ** (1.0 / max(self.geometry.dim, 1))) - avg_density = ideal_distance / max(mean_distance, BOUNDARY_FUZZ) - avg_force = float(np.linalg.norm(displacements[index])) + avg_density = float(force_summary.point_density[index]) + avg_force = float(force_summary.point_average_force[index]) fate = handler(self.rng, (avg_density, avg_force), thresh_add, thresh_del) if fate == PointFate.ADD_ANOTHER: @@ -126,6 +191,11 @@ def _attempt_add_delete_points(self, displacements: FloatArray) -> None: elif fate == PointFate.DELETE and len(self.points) - len(removals) > self.geometry.dim + 1: removals.append(index) + return additions, removals + + def _apply_point_changes(self, additions: list[FloatArray], removals: list[int]) -> None: + """Mutate the point cloud according to the evaluated add/remove plan.""" + if removals: keep = np.ones(len(self.points), dtype=bool) keep[np.asarray(removals, dtype=int)] = False @@ -136,69 +206,124 @@ def _attempt_add_delete_points(self, displacements: FloatArray) -> None: additions_arr = _dedupe_points(np.asarray(additions, dtype=float)) if len(additions_arr) > 0: self.points = np.vstack((self.points, additions_arr)) + addition_states = _classify_dynamic_states(self.geometry, additions_arr, self.a0) self.states = np.concatenate( ( self.states, - np.full(len(additions_arr), STATE_MOBILE, dtype=int), + addition_states, ) ) - def _step_once(self) -> None: - """Execute one relaxation step over all currently mobile points.""" - - if len(self.points) < self.geometry.dim + 1: - self.finished = True - return + def _compute_constrained_displacement( + self, + point: FloatArray, + state: int, + displacement: FloatArray, + max_norm: float, + ) -> FloatArray: + """Return a clipped, boundary-aware candidate position for one point.""" + + bounded_displacement = np.array(displacement, copy=True) + if state == STATE_BOUNDARY: + normal = self.geometry.boundary_normal(point) + normal_norm = float(np.linalg.norm(normal)) + if normal_norm > BOUNDARY_FUZZ: + bounded_displacement -= float(np.dot(bounded_displacement, normal)) * (normal / normal_norm) + + disp_norm = float(np.linalg.norm(bounded_displacement)) + if disp_norm > max_norm > 0.0: + bounded_displacement *= max_norm / disp_norm + + candidate = point + bounded_displacement + if state == STATE_BOUNDARY: + old_boundary_distance = self.geometry.boundary_distance(point) + new_boundary_distance = self.geometry.boundary_distance(candidate) + if new_boundary_distance > old_boundary_distance + self.boundary_drift_tolerance: + return np.asarray(point, dtype=float) + + if self.geometry.classify_points(candidate.reshape(1, -1))[0] < 0: + candidate = self.geometry.project_segment_to_domain(point, candidate) + + if self.geometry.classify_points(candidate.reshape(1, -1))[0] >= 0: + return np.asarray(candidate, dtype=float) + return np.asarray(point, dtype=float) + + def _apply_relaxation_displacements( + self, + force_summary: ForceSummary, + time_step: float, + ) -> float: + """Move dynamic points according to the current force field.""" - self._rebuild_mesh() - neighbors = self._neighbor_map() new_points = np.array(self.points, copy=True) - displacements = np.zeros_like(new_points) max_displacement = 0.0 - + max_norm = self.a0 * float(self.params.get("controller_movement_max_freedom", 3.0)) * 0.05 for index, state in enumerate(self.states): - if state != STATE_MOBILE: - continue - - neigh = neighbors[index] - if not neigh: + if state not in (STATE_MOBILE, STATE_BOUNDARY): continue point = self.points[index] - neigh_coords = self.points[np.asarray(neigh, dtype=int)] - deltas = neigh_coords - point - distances = np.linalg.norm(deltas, axis=1) - nonzero = distances > BOUNDARY_FUZZ - if not np.any(nonzero): - continue + displacement = np.array(force_summary.total_forces[index], copy=True) * time_step + candidate = self._compute_constrained_displacement(point, int(state), displacement, max_norm) + new_points[index] = candidate + max_displacement = max(max_displacement, float(np.linalg.norm(candidate - point))) + self.points = new_points + return max_displacement + + def _check_convergence(self) -> None: + """Update the engine state when a convergence condition is satisfied.""" + + if ( + self.step >= self.min_equilibrium_steps + and self.last_max_displacement <= self.tolerated_rel_move + ): + log.debug( + "Relaxation finished by movement convergence at step %d (rel_move=%g, effective_force=%g)", + self.step, + self.last_max_displacement, + self.last_max_effective_force, + ) + self.finished = True + return - deltas = deltas[nonzero] - distances = distances[nonzero] - density_here = self.geometry.density_at(point) - ideal_distance = self.a0 / (density_here ** (1.0 / max(self.geometry.dim, 1))) - scale = (distances - ideal_distance) / distances - displacement = np.mean(deltas * scale[:, None], axis=0) - displacement *= self.time_step_scale - max_norm = self.a0 * float(self.params.get("controller_movement_max_freedom", 3.0)) * 0.05 - disp_norm = float(np.linalg.norm(displacement)) - if disp_norm > max_norm > 0.0: - displacement *= max_norm / disp_norm - - candidate = point + displacement - if self.geometry.classify_points(candidate.reshape(1, -1))[0] >= 0: - new_points[index] = candidate - displacements[index] = displacement - max_displacement = max(max_displacement, float(np.linalg.norm(displacement))) + if ( + self.step >= self.min_equilibrium_steps + and self.last_max_effective_force <= BOUNDARY_FUZZ + ): + log.debug( + "Relaxation finished by effective-force equilibrium at step %d (effective_force=%g, rel_move=%g)", + self.step, + self.last_max_effective_force, + self.last_max_displacement, + ) + self.finished = True - self.points = new_points + def _step_once(self) -> None: + """Execute one relaxation step over all currently mobile or boundary points.""" + + if len(self.points) < self.geometry.dim + 1: + self.finished = True + return + + force_summary = compute_forces( + self.points, + self.states, + self.geometry, + self.a0, + self.params, + self.step, + ) + time_step = self._effective_time_step(force_summary) + max_displacement = self._apply_relaxation_displacements(force_summary, time_step) self.last_max_displacement = max_displacement / max(self.a0, BOUNDARY_FUZZ) + self.last_max_effective_force = force_summary.max_effective_force + self._refresh_boundary_states() if self.step > 0 and self.step % 5 == 0: - self._attempt_add_delete_points(displacements) + self._attempt_add_delete_points(force_summary) self._rebuild_mesh() - if self.last_max_displacement <= self.tolerated_rel_move: - self.finished = True + self._check_convergence() def callback_payload(self) -> list[list[Any]]: """Return the callback payload for the engine's current mesh snapshot.""" diff --git a/src/nmesh/mesher/relaxation/forces.py b/src/nmesh/mesher/relaxation/forces.py new file mode 100644 index 0000000..ea728b2 --- /dev/null +++ b/src/nmesh/mesher/relaxation/forces.py @@ -0,0 +1,777 @@ +"""Force calculations for the relaxation meshing pipeline.""" + +from __future__ import annotations + +import math +from dataclasses import dataclass +from itertools import combinations +from typing import Any + +import numpy as np + +from ..meshing_parameters import ( + default_boundary_node_force_fun, + default_initial_relaxation_weight, + default_relaxation_force_fun, +) +from ._constants import ( + BOUNDARY_FUZZ, + DENSITY_EPSILON, + STATE_BOUNDARY, + STATE_FIXED, + STATE_MOBILE, + STATE_SIMPLE, +) +from ._types import FloatArray +from .geometry import FemGeometry +from .jit import accumulate_neighbor_forces_default +from .topology import _simplex_measures, _triangulate_points + + +@dataclass(frozen=True, slots=True) +class ForceSummary: + """Container for the per-step force calculation results.""" + + total_forces: FloatArray + neighbor_map: list[list[int]] + simplices: np.ndarray + point_density: FloatArray + point_average_force: FloatArray + point_effective_force: FloatArray + max_effective_force: float + + +@dataclass(frozen=True, slots=True) +class ForceParameters: + """Resolved controller parameters for one relaxation-force evaluation.""" + + shape_force_scale: float + volume_force_scale: float + neigh_force_scale: float + irrel_force_scale: float + sliver_correction: float + relaxation_weight: float + force_fun: Any + boundary_force_fun: Any + + +def _corner_force_threshold(dim: int) -> float: + """Return the solid-angle threshold below which corner suppression applies.""" + + if dim == 2: + return 0.75 * math.pi + if dim == 3: + return 0.85 * math.pi + return -math.inf + + +def _build_neighbor_pairs(simplices: np.ndarray) -> np.ndarray: + """Return sorted unique point-pair edges induced by the simplices.""" + + if len(simplices) == 0: + return np.empty((0, 2), dtype=np.int64) + + pairs: set[tuple[int, int]] = set() + for simplex in simplices: + for left, right in combinations(simplex.tolist(), 2): + i = int(left) + j = int(right) + pairs.add((i, j) if i <= j else (j, i)) + return np.asarray(sorted(pairs), dtype=np.int64) + + +def _is_dynamic_state(state: int) -> bool: + """Return whether the state participates in relaxation movement.""" + + return state in (STATE_MOBILE, STATE_BOUNDARY) + + +def _is_boundary_interaction(state_a: int, state_b: int) -> bool: + """Return whether a neighbor interaction should use the boundary law.""" + + return ( + state_a in (STATE_FIXED, STATE_BOUNDARY, STATE_SIMPLE) + or state_b in (STATE_FIXED, STATE_BOUNDARY, STATE_SIMPLE) + ) + + +def _regular_simplex_volume(edge_length: float, dim: int) -> float: + """Return the volume of a regular simplex with the supplied edge length.""" + + if dim <= 0: + return edge_length + numerator = edge_length**dim * math.sqrt(dim + 1.0) + denominator = math.factorial(dim) * math.sqrt(2.0**dim) + return numerator / denominator + + +def _sphere_volume(radius: float, dim: int) -> float: + """Return the d-dimensional volume of a sphere.""" + + return (math.pi ** (0.5 * dim) / math.gamma(0.5 * dim + 1.0)) * (radius**dim) + + +def _neighbor_forces_python( + points: FloatArray, + states: np.ndarray, + point_densities: FloatArray, + neighbor_map: list[list[int]], + dim: int, + a0: float, + neigh_force_scale: float, + force_fun: Any, + boundary_force_fun: Any, +) -> tuple[FloatArray, FloatArray, np.ndarray]: + """Fallback neighbor-force path for custom Python force callbacks.""" + + total_forces = np.zeros((len(points), dim), dtype=float) + neighbor_force_sums = np.zeros(len(points), dtype=float) + neighbor_force_counts = np.zeros(len(points), dtype=int) + + for left, neighbors in enumerate(neighbor_map): + point_left = points[left] + state_left = int(states[left]) + for right in neighbors: + if right <= left: + continue + contribution, scalar_force = _apply_pairwise_force( + point_left=point_left, + point_right=points[right], + state_left=state_left, + state_right=int(states[right]), + density_left=float(point_densities[left]), + density_right=float(point_densities[right]), + dim=dim, + a0=a0, + neigh_force_scale=neigh_force_scale, + force_fun=force_fun, + boundary_force_fun=boundary_force_fun, + ) + if scalar_force is None: + continue + neighbor_force_sums[left] += abs(scalar_force) + neighbor_force_sums[right] += abs(scalar_force) + neighbor_force_counts[left] += 1 + neighbor_force_counts[right] += 1 + if contribution is None: + continue + total_forces[left] += contribution + total_forces[right] -= contribution + + return total_forces, neighbor_force_sums, neighbor_force_counts + + +def _apply_pairwise_force( + *, + point_left: FloatArray, + point_right: FloatArray, + state_left: int, + state_right: int, + density_left: float, + density_right: float, + dim: int, + a0: float, + neigh_force_scale: float, + force_fun: Any, + boundary_force_fun: Any, +) -> tuple[FloatArray | None, float | None]: + """Compute the pairwise neighbor contribution for one point pair.""" + + if not _is_dynamic_state(state_left) and not _is_dynamic_state(state_right): + return None, None + + delta = point_right - point_left + true_distance = float(np.linalg.norm(delta)) + if true_distance <= DENSITY_EPSILON: + return None, None + + avg_density = 0.5 * (density_left + density_right) + inv_length_scale = (avg_density ** (1.0 / max(dim, 1))) / max(a0, DENSITY_EPSILON) + reduced_distance = true_distance * inv_length_scale + scalar_force = float( + (boundary_force_fun if _is_boundary_interaction(state_left, state_right) else force_fun)( + reduced_distance + ) + ) + scaled_force = neigh_force_scale * abs(scalar_force) + if scalar_force == 0.0: + return None, scaled_force + return neigh_force_scale * (-scalar_force) * delta, scaled_force + + +def _extract_force_parameters(params: dict[str, Any], step: int) -> ForceParameters: + """Resolve the mesher controller parameters used by ``compute_forces``.""" + + settling_steps = int(params.get("controller_initial_settling_steps", 100)) + relaxation_weight_fun = params.get("initial_relaxation_weight_fun", default_initial_relaxation_weight) + return ForceParameters( + shape_force_scale=float(params.get("controller_shape_force_scale", 0.1)), + volume_force_scale=float(params.get("controller_volume_force_scale", 0.0)), + neigh_force_scale=float(params.get("controller_neigh_force_scale", 1.0)), + irrel_force_scale=float(params.get("controller_irrel_elem_force_scale", 1.0)), + sliver_correction=float(params.get("controller_sliver_correction", 1.0)), + relaxation_weight=float(relaxation_weight_fun(step, settling_steps, 0.0, 1.0)), + force_fun=params.get("relaxation_force_fun", default_relaxation_force_fun), + boundary_force_fun=params.get("boundary_node_force_fun", default_boundary_node_force_fun), + ) + + +def _build_neighbor_map(point_count: int, simplices: np.ndarray) -> list[list[int]]: + """Build undirected point adjacency from simplices.""" + + neighbors: list[set[int]] = [set() for _ in range(point_count)] + for simplex in simplices: + for left, right in combinations(simplex.tolist(), 2): + i = int(left) + j = int(right) + neighbors[i].add(j) + neighbors[j].add(i) + return [sorted(group) for group in neighbors] + + +def _compute_neighbor_forces( + points: FloatArray, + states: np.ndarray, + simplices: np.ndarray, + neighbor_map: list[list[int]], + point_densities: FloatArray, + dim: int, + a0: float, + config: ForceParameters, +) -> tuple[FloatArray, FloatArray, np.ndarray]: + """Compute neighbor-force contributions through the JIT or Python path.""" + + point_count = len(points) + if len(simplices) == 0: + return ( + np.zeros((point_count, dim), dtype=float), + np.zeros(point_count, dtype=float), + np.zeros(point_count, dtype=int), + ) + + if ( + config.force_fun is default_relaxation_force_fun + and config.boundary_force_fun is default_boundary_node_force_fun + ): + neighbor_pairs = _build_neighbor_pairs(simplices) + return accumulate_neighbor_forces_default( + np.asarray(points, dtype=np.float64), + np.asarray(states, dtype=np.int64), + point_densities.astype(np.float64, copy=False), + neighbor_pairs, + float(a0), + config.neigh_force_scale, + ) + + return _neighbor_forces_python( + points, + states, + point_densities, + neighbor_map, + dim, + a0, + config.neigh_force_scale, + config.force_fun, + config.boundary_force_fun, + ) + + +def _vertex_angle(points: FloatArray, simplex: np.ndarray, local_index: int, dim: int) -> float: + """Return the angle or solid angle covered by a simplex at one vertex.""" + + vertex = points[int(simplex[local_index])] + others = [ + points[int(simplex[index])] - vertex + for index in range(len(simplex)) + if index != local_index + ] + + if dim == 1: + return math.pi + + if dim == 2: + first = others[0] + second = others[1] + denom = float(np.linalg.norm(first) * np.linalg.norm(second)) + if denom <= DENSITY_EPSILON: + return 0.0 + cosine = float(np.clip(np.dot(first, second) / denom, -1.0, 1.0)) + return float(math.acos(cosine)) + + if dim == 3: + a, b, c = others + numer = abs(float(np.dot(a, np.cross(b, c)))) + denom = ( + float(np.linalg.norm(a) * np.linalg.norm(b) * np.linalg.norm(c)) + + float(np.dot(a, b) * np.linalg.norm(c)) + + float(np.dot(a, c) * np.linalg.norm(b)) + + float(np.dot(b, c) * np.linalg.norm(a)) + ) + if numer <= DENSITY_EPSILON and denom <= DENSITY_EPSILON: + return 0.0 + return float(2.0 * math.atan2(numer, max(denom, DENSITY_EPSILON))) + + return 1.0 + + +def _shape_force_matrix(vertices: FloatArray, dim: int) -> FloatArray: + """Build the covariance-derived shape force matrix for one simplex.""" + + covariance = np.asarray(vertices.T @ vertices, dtype=float) + covariance /= max(float(len(vertices) - 1), 1.0) + + eigenvalues, eigenvectors = np.linalg.eigh(covariance) + eigenvectors = np.asarray(eigenvectors, dtype=float) + safe_eigenvalues = np.clip(eigenvalues.astype(float, copy=False), DENSITY_EPSILON, None) + product = float(np.prod(safe_eigenvalues)) + if product <= DENSITY_EPSILON: + return np.zeros((dim, dim), dtype=float) + + scaled = safe_eigenvalues * ((1.0 / product) ** (1.0 / max(dim, 1))) + transformed = np.asarray(np.diag(-np.log(np.clip(scaled, DENSITY_EPSILON, None))), dtype=float) + return np.asarray(eigenvectors @ transformed @ eigenvectors.T, dtype=float) + + +def _project_force_to_tangent( + geometry: FemGeometry, + point: FloatArray, + force: FloatArray, +) -> FloatArray: + """Remove the normal component of a force at a boundary point.""" + + normal = geometry.boundary_normal(point) + normal_norm = float(np.linalg.norm(normal)) + if normal_norm <= BOUNDARY_FUZZ: + return np.asarray(force, dtype=float) + unit_normal = normal / normal_norm + return np.asarray(force, dtype=float) - float(np.dot(force, unit_normal)) * unit_normal + + +def _simplex_force_contribution( + simplex_points: FloatArray, + dim: int, + density_here: float, + a0: float, + relaxation_weight: float, + shape_force_scale: float, + volume_force_scale: float, + sliver_correction: float, + volume: float, +) -> FloatArray: + """Compute the shape and volume force contribution for one simplex.""" + + center = np.mean(simplex_points, axis=0) + vertices = simplex_points - center + ideal_edge_length = a0 / (density_here ** (1.0 / max(dim, 1))) + ideal_volume = _regular_simplex_volume(ideal_edge_length, dim) + return _compute_volume_forces( + vertices, + relaxation_weight, + volume_force_scale, + volume, + ideal_volume, + ) + _compute_shape_forces( + vertices=vertices, + dim=dim, + relaxation_weight=relaxation_weight, + shape_force_scale=shape_force_scale, + sliver_correction=sliver_correction, + volume=volume, + ideal_volume=ideal_volume, + ) + + +def _compute_volume_forces( + vertices: FloatArray, + relaxation_weight: float, + volume_force_scale: float, + volume: float, + ideal_volume: float, +) -> FloatArray: + """Return the isotropic volume-restoring forces for one simplex.""" + + forces = np.zeros_like(vertices) + if volume_force_scale <= 0.0 or volume <= DENSITY_EPSILON or ideal_volume <= DENSITY_EPSILON: + return forces + volume_factor = volume_force_scale * math.log(ideal_volume / volume) + return forces + relaxation_weight * volume_factor * vertices + + +def _compute_shape_forces( + *, + vertices: FloatArray, + dim: int, + relaxation_weight: float, + shape_force_scale: float, + sliver_correction: float, + volume: float, + ideal_volume: float, +) -> FloatArray: + """Return the covariance-derived shape-correction forces for one simplex.""" + + forces = np.zeros_like(vertices) + if shape_force_scale <= 0.0: + return forces + + shape_matrix = _shape_force_matrix(vertices, dim) + if not np.any(shape_matrix): + return forces + + vol_correction = _shape_volume_correction(volume, ideal_volume, dim) + for index, vertex in enumerate(vertices): + norm = float(np.linalg.norm(vertex)) + if norm <= DENSITY_EPSILON: + continue + raw_force = shape_matrix @ vertex + normal = vertex / norm + projection = float(np.dot(raw_force, normal)) + angular = raw_force - projection * normal + forces[index] += relaxation_weight * vol_correction * shape_force_scale * angular + if projection > 0.0 and sliver_correction > 0.0: + longitudinal = raw_force - angular + forces[index] += ( + shape_force_scale + * sliver_correction + * relaxation_weight + * max(vol_correction - 1.0, 0.0) + * longitudinal + ) + return forces + + +def _shape_volume_correction(volume: float, ideal_volume: float, dim: int) -> float: + """Return the volume-dependent multiplier applied to shape forces.""" + + if volume <= DENSITY_EPSILON or ideal_volume <= DENSITY_EPSILON: + return 0.0 + offset = 0.0 if dim <= 2 else 1.0 + return max(1.0, offset + math.log(ideal_volume / volume)) + + +def _simplex_incidence_data( + points: FloatArray, + simplices: np.ndarray, + dim: int, +) -> tuple[FloatArray, list[list[int]], FloatArray]: + """Collect simplex measures plus per-point incident simplex and angle data.""" + + simplex_measures = _simplex_measures(points, simplices, dim) + incident_simplices: list[list[int]] = [[] for _ in range(len(points))] + angle_sums = np.zeros(len(points), dtype=float) + for simplex_index, simplex in enumerate(simplices): + for local_index, point_index in enumerate(simplex.tolist()): + point_id = int(point_index) + incident_simplices[point_id].append(simplex_index) + angle_sums[point_id] += _vertex_angle(points, simplex, local_index, dim) + return simplex_measures, incident_simplices, angle_sums + + +def _compute_simplex_forces( + total_forces: FloatArray, + points: FloatArray, + states: np.ndarray, + geometry: FemGeometry, + simplices: np.ndarray, + point_densities: FloatArray, + simplex_measures: FloatArray, + angle_sums: FloatArray, + dim: int, + a0: float, + config: ForceParameters, +) -> None: + """Accumulate simplex shape, volume, and irrelevant-element forces.""" + + if len(simplices) == 0: + return + if ( + config.shape_force_scale <= 0.0 + and config.volume_force_scale <= 0.0 + and config.irrel_force_scale <= 0.0 + ): + return + + centroid_regions = geometry.classify_points(np.mean(points[simplices], axis=1)) + suppress_corner_forces = angle_sums < _corner_force_threshold(dim) + for simplex_index, simplex in enumerate(simplices): + simplex_points = points[simplex] + simplex_force = np.zeros_like(simplex_points) + volume = float(simplex_measures[simplex_index]) + if centroid_regions[simplex_index] >= 0 and volume > BOUNDARY_FUZZ: + density_here = float(np.mean(point_densities[np.asarray(simplex, dtype=int)])) + simplex_force = _simplex_force_contribution( + simplex_points=simplex_points, + dim=dim, + density_here=density_here, + a0=a0, + relaxation_weight=config.relaxation_weight, + shape_force_scale=config.shape_force_scale, + volume_force_scale=config.volume_force_scale, + sliver_correction=config.sliver_correction, + volume=volume, + ) + elif config.irrel_force_scale > 0.0: + center = np.mean(simplex_points, axis=0) + for local_index, point_index in enumerate(simplex.tolist()): + if not _is_dynamic_state(int(states[point_index])): + continue + simplex_force[local_index] = config.irrel_force_scale * (center - simplex_points[local_index]) + + for local_index, point_index in enumerate(simplex.tolist()): + point_id = int(point_index) + if suppress_corner_forces[point_id]: + continue + total_forces[point_id] += simplex_force[local_index] + + +def _finalize_force_summary( + total_forces: FloatArray, + neighbor_map: list[list[int]], + simplices: np.ndarray, + point_densities: FloatArray, + neighbor_force_sums: FloatArray, + neighbor_force_counts: np.ndarray, + simplex_measures: FloatArray, + incident_simplices: list[list[int]], + angle_sums: FloatArray, + points: FloatArray, + states: np.ndarray, + geometry: FemGeometry, + a0: float, + dim: int, +) -> ForceSummary: + """Assemble the per-point density and effective-force metrics.""" + + point_count = len(points) + point_density, point_average_force, point_effective_force = _empty_force_arrays(point_count) + + if point_count == 0: + return _make_force_summary( + total_forces, + neighbor_map, + simplices, + point_density, + point_average_force, + point_effective_force, + ) + + boundary_mask = geometry.boundary_mask(points, tolerance=max(0.05 * a0, BOUNDARY_FUZZ)) + full_angle = _full_angle(dim) + + for point_index, point in enumerate(points): + point_average_force[point_index] = _average_neighbor_force( + neighbor_force_sums[point_index], + int(neighbor_force_counts[point_index]), + ) + corrected_volume = _corrected_point_volume( + simplex_measures=simplex_measures, + incident=incident_simplices[point_index], + angle=float(angle_sums[point_index]), + dim=dim, + a0=a0, + is_boundary=bool(boundary_mask[point_index]), + full_angle=full_angle, + ) + density_here = float(point_densities[point_index]) + point_density[point_index] = _point_density_ratio( + density_here=density_here, + corrected_volume=corrected_volume, + a0=a0, + dim=dim, + ) + point_effective_force[point_index] = _point_effective_force( + total_force=total_forces[point_index], + state=int(states[point_index]), + geometry=geometry, + point=point, + is_boundary=bool(boundary_mask[point_index]), + density_here=density_here, + a0=a0, + dim=dim, + ) + + return _make_force_summary( + total_forces, + neighbor_map, + simplices, + point_density, + point_average_force, + point_effective_force, + ) + + +def _empty_force_arrays(point_count: int) -> tuple[FloatArray, FloatArray, FloatArray]: + """Allocate the per-point summary arrays used by ``_finalize_force_summary``.""" + + return ( + np.ones(point_count, dtype=float), + np.zeros(point_count, dtype=float), + np.zeros(point_count, dtype=float), + ) + + +def _make_force_summary( + total_forces: FloatArray, + neighbor_map: list[list[int]], + simplices: np.ndarray, + point_density: FloatArray, + point_average_force: FloatArray, + point_effective_force: FloatArray, +) -> ForceSummary: + """Construct a ``ForceSummary`` from the completed per-point arrays.""" + + return ForceSummary( + total_forces=total_forces, + neighbor_map=neighbor_map, + simplices=simplices, + point_density=point_density, + point_average_force=point_average_force, + point_effective_force=point_effective_force, + max_effective_force=float(np.max(point_effective_force, initial=0.0)), + ) + + +def _full_angle(dim: int) -> float: + """Return the complete angular measure used for Voronoi-style correction.""" + + if dim == 1: + return 1.5 + if dim == 2: + return 2.0 * math.pi + if dim == 3: + return 4.0 * math.pi + return 1.0 + + +def _average_neighbor_force(force_sum: float, force_count: int) -> float: + """Return the mean absolute neighbor-force magnitude for a point.""" + + if force_count <= 0: + return 0.0 + return force_sum / float(force_count) + + +def _corrected_point_volume( + *, + simplex_measures: FloatArray, + incident: list[int], + angle: float, + dim: int, + a0: float, + is_boundary: bool, + full_angle: float, +) -> float: + """Return the local corrected control volume for one point.""" + + if incident: + local_volume = float(np.sum(simplex_measures[np.asarray(incident, dtype=int)])) / float(dim + 1) + else: + local_volume = _sphere_volume(0.5 * a0, dim) + + if angle <= DENSITY_EPSILON: + return 1.0e-4 * _sphere_volume(0.5 * a0, dim) + + correction = ( + full_angle / angle + if dim in (1, 2, 3) + else max(1.0, float(dim + 1) / float(len(incident) or 1)) + ) + corrected_volume = local_volume * correction + if is_boundary: + corrected_volume *= 1.2 + return corrected_volume + + +def _point_density_ratio( + *, + density_here: float, + corrected_volume: float, + a0: float, + dim: int, +) -> float: + """Return the desired-to-actual local density ratio for one point.""" + + effective_rod_length = a0 / (density_here ** (1.0 / max(dim, 1))) + ideal_local_volume = _sphere_volume(0.5 * effective_rod_length, dim) + return ideal_local_volume / max(corrected_volume, DENSITY_EPSILON) + + +def _point_effective_force( + *, + total_force: FloatArray, + state: int, + geometry: FemGeometry, + point: FloatArray, + is_boundary: bool, + density_here: float, + a0: float, + dim: int, +) -> float: + """Return the normalized effective force magnitude for one point.""" + + effective_force = np.asarray(total_force, dtype=float) + if state == STATE_BOUNDARY or is_boundary: + effective_force = _project_force_to_tangent(geometry, point, effective_force) + return ( + float(np.linalg.norm(effective_force)) + * (density_here ** (1.0 / max(dim, 1))) + / max(a0, DENSITY_EPSILON) + ) + + +def compute_forces( + points: FloatArray, + states: np.ndarray, + geometry: FemGeometry, + a0: float, + params: dict[str, Any], + step: int, +) -> ForceSummary: + """Compute neighbor, shape, volume, and irrelevant-element forces.""" + + point_count = len(points) + dim = geometry.dim + simplices = _triangulate_points(points, dim, states) + neighbor_map = _build_neighbor_map(point_count, simplices) + config = _extract_force_parameters(params, step) + point_densities = np.asarray([geometry.density_at(point) for point in points], dtype=float) + + total_forces, neighbor_force_sums, neighbor_force_counts = _compute_neighbor_forces( + points, + states, + simplices, + neighbor_map, + point_densities, + dim, + a0, + config, + ) + simplex_measures, incident_simplices, angle_sums = _simplex_incidence_data(points, simplices, dim) + _compute_simplex_forces( + total_forces, + points, + states, + geometry, + simplices, + point_densities, + simplex_measures, + angle_sums, + dim, + a0, + config, + ) + return _finalize_force_summary( + total_forces=total_forces, + neighbor_map=neighbor_map, + simplices=simplices, + point_densities=point_densities, + neighbor_force_sums=neighbor_force_sums, + neighbor_force_counts=neighbor_force_counts, + simplex_measures=simplex_measures, + incident_simplices=incident_simplices, + angle_sums=angle_sums, + points=points, + states=states, + geometry=geometry, + a0=a0, + dim=dim, + ) diff --git a/src/nmesh/mesher/relaxation/geometry.py b/src/nmesh/mesher/relaxation/geometry.py index d039ee3..da99c21 100644 --- a/src/nmesh/mesher/relaxation/geometry.py +++ b/src/nmesh/mesher/relaxation/geometry.py @@ -2,6 +2,7 @@ from __future__ import annotations +import math from dataclasses import dataclass from typing import Any @@ -24,6 +25,7 @@ class FemGeometry: density_fun: DensityFunction region_ids: tuple[int, ...] region_functions: tuple[RegionFunction, ...] + region_bodies: tuple[Body | None, ...] piece_hints: tuple[FloatArray, ...] mesh_exterior: bool @@ -57,6 +59,116 @@ def density_at(self, point: FloatArray) -> float: return max(float(self.density_fun(point)), DENSITY_EPSILON) + def boundary_mask(self, points: FloatArray, tolerance: float) -> np.ndarray: + """Return a mask for points that lie close to a region or box boundary.""" + + if points.size == 0: + return np.empty(0, dtype=bool) + + coords = np.asarray(points, dtype=float) + mask = np.zeros(len(coords), dtype=bool) + for body in self.region_bodies: + if body is None: + continue + mask |= np.abs(np.asarray(body.evaluate(coords), dtype=float)) <= tolerance + + if self.mesh_exterior: + mask |= np.any( + np.isclose(coords, self.bbox_min, atol=tolerance, rtol=0.0) + | np.isclose(coords, self.bbox_max, atol=tolerance, rtol=0.0), + axis=1, + ) + + return mask & self.points_in_bbox(coords) + + def boundary_distance(self, point: FloatArray) -> float: + """Return the smallest available distance-like scalar to a boundary.""" + + coord = np.asarray(point, dtype=float) + distances: list[float] = [] + for body in self.region_bodies: + if body is None: + continue + distances.append(abs(float(body.evaluate(coord)))) + + if self.mesh_exterior: + distances.extend( + min(abs(float(coord[axis] - self.bbox_min[axis])), abs(float(coord[axis] - self.bbox_max[axis]))) + for axis in range(self.dim) + ) + + if not distances: + return math.inf + return min(distances) + + def boundary_normal(self, point: FloatArray) -> FloatArray: + """Estimate a unit normal for the closest relevant boundary.""" + + coord = np.asarray(point, dtype=float) + best_distance = math.inf + best_normal = np.zeros(self.dim, dtype=float) + extent = max(float(np.max(self.bbox_max - self.bbox_min)), 1.0) + epsilon = max(BOUNDARY_FUZZ * 10.0, extent * 1.0e-6) + + for body in self.region_bodies: + if body is None: + continue + distance = abs(float(body.evaluate(coord))) + if distance > best_distance: + continue + gradient = np.zeros(self.dim, dtype=float) + for axis in range(self.dim): + offset = np.zeros(self.dim, dtype=float) + offset[axis] = epsilon + forward = float(body.evaluate(coord + offset)) + backward = float(body.evaluate(coord - offset)) + gradient[axis] = (forward - backward) / (2.0 * epsilon) + norm = float(np.linalg.norm(gradient)) + if norm <= DENSITY_EPSILON: + continue + best_distance = distance + best_normal = gradient / norm + + if self.mesh_exterior: + for axis in range(self.dim): + dist_min = abs(float(coord[axis] - self.bbox_min[axis])) + if dist_min < best_distance: + best_distance = dist_min + best_normal = np.zeros(self.dim, dtype=float) + best_normal[axis] = 1.0 + dist_max = abs(float(coord[axis] - self.bbox_max[axis])) + if dist_max < best_distance: + best_distance = dist_max + best_normal = np.zeros(self.dim, dtype=float) + best_normal[axis] = -1.0 + + return best_normal + + def project_segment_to_domain( + self, + start: FloatArray, + end: FloatArray, + *, + iterations: int = 18, + ) -> FloatArray: + """Project an outside point back into the domain along a segment.""" + + lower = np.asarray(start, dtype=float).copy() + upper = np.asarray(end, dtype=float).copy() + if self.classify_points(lower.reshape(1, -1))[0] < 0: + return lower + if self.classify_points(upper.reshape(1, -1))[0] >= 0: + return upper + + for _ in range(iterations): + middle = 0.5 * (lower + upper) + if self.classify_points(middle.reshape(1, -1))[0] >= 0: + lower = middle + else: + upper = middle + + return lower + def _raw_mesh_points(raw_mesh: RawMesh, dim: int) -> FloatArray: """Return mesh points as a validated floating-point array.""" @@ -88,24 +200,17 @@ def _dedupe_hint_points(points: FloatArray) -> FloatArray: return points[np.asarray(keep_indices, dtype=int)] -def fem_geometry_from_bodies( - bounding_box: tuple[FloatArray, FloatArray], +def _process_bodies_and_hints( bodies: list[Body], hints: list[list[Any]], - *, - density: str | DensityFunction | None = None, - mesh_exterior: bool = False, -) -> FemGeometry: - """Construct the geometry bundle consumed by the Python meshing engine.""" - - bbox_min = np.asarray(bounding_box[0], dtype=float) - bbox_max = np.asarray(bounding_box[1], dtype=float) - dim = int(len(bbox_min)) - density_fun = _compile_density_function(density) + dim: int, +) -> tuple[list[RegionFunction], list[FloatArray], list[int], list[Body | None]]: + """Convert bodies and hint meshes into region predicates and hint blocks.""" body_functions: list[RegionFunction] = [] piece_hints: list[FloatArray] = [] region_ids: list[int] = [] + region_bodies: list[Body | None] = [] next_region_id = 1 for body in bodies: @@ -115,6 +220,7 @@ def fem_geometry_from_bodies( ) piece_hints.append(np.empty((0, dim), dtype=float)) region_ids.append(next_region_id) + region_bodies.append(body) next_region_id += 1 for hint_mesh, hint_body in hints: @@ -128,10 +234,26 @@ def fem_geometry_from_bodies( ) piece_hints.append(hint_points) region_ids.append(next_region_id) + region_bodies.append(hint_body) next_region_id += 1 + return body_functions, piece_hints, region_ids, region_bodies + + +def _build_region_lists( + bbox_min: FloatArray, + bbox_max: FloatArray, + body_functions: list[RegionFunction], + region_ids: list[int], + region_bodies: list[Body | None], + *, + mesh_exterior: bool, +) -> tuple[list[RegionFunction], list[int], list[Body | None]]: + """Build the ordered region lists, including the exterior region when requested.""" + region_functions: list[RegionFunction] = [] ordered_region_ids: list[int] = [] + ordered_region_bodies: list[Body | None] = [] if mesh_exterior: @@ -150,9 +272,41 @@ def exterior(points: FloatArray) -> np.ndarray: region_functions.append(exterior) ordered_region_ids.append(0) + ordered_region_bodies.append(None) region_functions.extend(body_functions) ordered_region_ids.extend(region_ids) + ordered_region_bodies.extend(region_bodies) + return region_functions, ordered_region_ids, ordered_region_bodies + + +def fem_geometry_from_bodies( + bounding_box: tuple[FloatArray, FloatArray], + bodies: list[Body], + hints: list[list[Any]], + *, + density: str | DensityFunction | None = None, + mesh_exterior: bool = False, +) -> FemGeometry: + """Construct the geometry bundle consumed by the Python meshing engine.""" + + bbox_min = np.asarray(bounding_box[0], dtype=float) + bbox_max = np.asarray(bounding_box[1], dtype=float) + dim = int(len(bbox_min)) + density_fun = _compile_density_function(density) + body_functions, piece_hints, region_ids, region_bodies = _process_bodies_and_hints( + bodies, + hints, + dim, + ) + region_functions, ordered_region_ids, ordered_region_bodies = _build_region_lists( + bbox_min, + bbox_max, + body_functions, + region_ids, + region_bodies, + mesh_exterior=mesh_exterior, + ) return FemGeometry( dim=dim, @@ -161,6 +315,7 @@ def exterior(points: FloatArray) -> np.ndarray: density_fun=density_fun, region_ids=tuple(ordered_region_ids), region_functions=tuple(region_functions), + region_bodies=tuple(ordered_region_bodies), piece_hints=tuple(piece_hints), mesh_exterior=bool(mesh_exterior), ) diff --git a/src/nmesh/mesher/relaxation/jit.py b/src/nmesh/mesher/relaxation/jit.py new file mode 100644 index 0000000..57748ea --- /dev/null +++ b/src/nmesh/mesher/relaxation/jit.py @@ -0,0 +1,102 @@ +"""Numba-accelerated kernels for the relaxation meshing pipeline.""" + +from __future__ import annotations + +import math + +import numpy as np +from numba import njit + +from ...utils.constants import MIN_DIVISION_MAGNITUDE +from ._constants import DENSITY_EPSILON, STATE_BOUNDARY, STATE_FIXED, STATE_MOBILE, STATE_SIMPLE + + +@njit(cache=True) +def _default_relaxation_force(reduced_distance: float) -> float: + """Return the legacy default mobile-mobile force law.""" + + if reduced_distance > 1.0: + return 0.0 + return 1.0 - reduced_distance + + +@njit(cache=True) +def _default_boundary_force(reduced_distance: float) -> float: + """Return the legacy default boundary interaction force law.""" + + if reduced_distance > 1.0: + return 0.0 + if reduced_distance < MIN_DIVISION_MAGNITUDE: + return 1.0e12 + return 1.0 / reduced_distance - 1.0 + + +@njit(cache=True) +def accumulate_neighbor_forces_default( + points: np.ndarray, + states: np.ndarray, + point_densities: np.ndarray, + neighbor_pairs: np.ndarray, + a0: float, + neigh_force_scale: float, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Accumulate default neighbor forces and controller statistics.""" + + point_count, dim = points.shape + total_forces = np.zeros((point_count, dim), dtype=np.float64) + neighbor_force_sums = np.zeros(point_count, dtype=np.float64) + neighbor_force_counts = np.zeros(point_count, dtype=np.int64) + inv_dim = 1.0 / float(max(dim, 1)) + a0_safe = max(a0, DENSITY_EPSILON) + + for pair_index in range(len(neighbor_pairs)): + left = int(neighbor_pairs[pair_index, 0]) + right = int(neighbor_pairs[pair_index, 1]) + state_left = int(states[left]) + state_right = int(states[right]) + left_dynamic = state_left == STATE_MOBILE or state_left == STATE_BOUNDARY + right_dynamic = state_right == STATE_MOBILE or state_right == STATE_BOUNDARY + if not left_dynamic and not right_dynamic: + continue + + true_distance_sq = 0.0 + for axis in range(dim): + delta = points[right, axis] - points[left, axis] + true_distance_sq += delta * delta + true_distance = math.sqrt(true_distance_sq) + if true_distance <= DENSITY_EPSILON: + continue + + avg_density = 0.5 * (point_densities[left] + point_densities[right]) + inv_length_scale = math.pow(avg_density, inv_dim) / a0_safe + reduced_distance = true_distance * inv_length_scale + boundary_interaction = ( + state_left == STATE_FIXED + or state_left == STATE_BOUNDARY + or state_left == STATE_SIMPLE + or state_right == STATE_FIXED + or state_right == STATE_BOUNDARY + or state_right == STATE_SIMPLE + ) + scalar_force = ( + _default_boundary_force(reduced_distance) + if boundary_interaction + else _default_relaxation_force(reduced_distance) + ) + + scaled_force = neigh_force_scale * abs(scalar_force) + neighbor_force_sums[left] += scaled_force + neighbor_force_sums[right] += scaled_force + neighbor_force_counts[left] += 1 + neighbor_force_counts[right] += 1 + + if scalar_force == 0.0: + continue + + force_factor = neigh_force_scale * (-scalar_force) + for axis in range(dim): + contribution = force_factor * (points[right, axis] - points[left, axis]) + total_forces[left, axis] += contribution + total_forces[right, axis] -= contribution + + return total_forces, neighbor_force_sums, neighbor_force_counts diff --git a/src/nmesh/mesher/relaxation/seeding.py b/src/nmesh/mesher/relaxation/seeding.py index 17ee460..0c5a5b1 100644 --- a/src/nmesh/mesher/relaxation/seeding.py +++ b/src/nmesh/mesher/relaxation/seeding.py @@ -6,7 +6,15 @@ import numpy as np -from ._constants import BOUNDARY_FUZZ, DEFAULT_RNG_SEED, DENSITY_EPSILON, STATE_FIXED, STATE_MOBILE, STATE_SIMPLE +from ._constants import ( + BOUNDARY_FUZZ, + DEFAULT_RNG_SEED, + DENSITY_EPSILON, + STATE_BOUNDARY, + STATE_FIXED, + STATE_MOBILE, + STATE_SIMPLE, +) from ._types import FloatArray from .geometry import FemGeometry @@ -120,70 +128,125 @@ def _keep_point_for_density(point: FloatArray, density_here: float, rng: np.rand return bool(rng.random() <= density_here) -def _prepare_initial_points( +def _classify_dynamic_states(geometry: FemGeometry, points: FloatArray, a0: float) -> np.ndarray: + """Return mobile or boundary states for the supplied movable points.""" + + if len(points) == 0: + return np.empty(0, dtype=int) + mask = geometry.boundary_mask(points, tolerance=max(0.05 * a0, BOUNDARY_FUZZ)) + states = np.full(len(points), STATE_MOBILE, dtype=int) + states[mask] = STATE_BOUNDARY + return states + + +def _filter_relevant_points(geometry: FemGeometry, points: FloatArray) -> FloatArray: + """Keep only points that belong to one of the meshed regions.""" + + if len(points) == 0: + return points + mask = geometry.classify_points(points) >= 0 + return points[mask] + + +def _select_generated_points( geometry: FemGeometry, a0: float, fixed_points: FloatArray, mobile_points: FloatArray, simply_points: FloatArray, - periodic: list[float] | list[bool], rng: np.random.Generator, -) -> tuple[FloatArray, np.ndarray]: - """Prepare the initial point cloud and point-state array for relaxation.""" +) -> FloatArray: + """Generate and thin candidate points that are not already seeded.""" - fixed_points, mobile_points, simply_points = _dedupe_fixed_mobile( - fixed_points, mobile_points, simply_points - ) + candidates = _grid_candidate_points(geometry, a0) + region_ids = geometry.classify_points(candidates) + candidates = candidates[region_ids >= 0] + if len(candidates) == 0: + return np.empty((0, geometry.dim), dtype=float) - def filter_relevant(points: FloatArray) -> FloatArray: - """Keep only points that belong to a meshed region.""" + candidate_keys = {_point_key(point) for point in fixed_points} + candidate_keys.update(_point_key(point) for point in mobile_points) + candidate_keys.update(_point_key(point) for point in simply_points) - if len(points) == 0: - return points - mask = geometry.classify_points(points) >= 0 - return points[mask] + selected_candidates: list[FloatArray] = [] + for point in candidates: + key = _point_key(point) + if key in candidate_keys: + continue + density_here = geometry.density_at(point) + if not _keep_point_for_density(point, density_here, rng): + continue + selected_candidates.append(point) + candidate_keys.add(key) - fixed_points = filter_relevant(fixed_points) - mobile_points = filter_relevant(mobile_points) - simply_points = filter_relevant(simply_points) + if not selected_candidates: + return np.empty((0, geometry.dim), dtype=float) + return np.asarray(selected_candidates, dtype=float) - candidates = _grid_candidate_points(geometry, a0) - region_ids = geometry.classify_points(candidates) - candidates = candidates[region_ids >= 0] - if len(candidates) > 0: - candidate_keys = {_point_key(point) for point in fixed_points} - candidate_keys.update(_point_key(point) for point in mobile_points) - candidate_keys.update(_point_key(point) for point in simply_points) +def _collect_hint_points(geometry: FemGeometry) -> FloatArray: + """Merge, deduplicate, and filter hint points from all geometry pieces.""" - selected_candidates: list[FloatArray] = [] - for point in candidates: - key = _point_key(point) - if key in candidate_keys: - continue - density_here = geometry.density_at(point) - if not _keep_point_for_density(point, density_here, rng): - continue - selected_candidates.append(point) - candidate_keys.add(key) - generated_points = ( - np.asarray(selected_candidates, dtype=float) - if selected_candidates - else np.empty((0, geometry.dim), dtype=float) + hint_points = [hint_points for hint_points in geometry.piece_hints if len(hint_points) > 0] + if not hint_points: + return np.empty((0, geometry.dim), dtype=float) + hint_block = _dedupe_points(np.vstack(hint_points)) + return _filter_relevant_points(geometry, hint_block) + + +def _apply_periodic_fixed_states( + states: np.ndarray, + all_points: FloatArray, + geometry: FemGeometry, + periodic: list[float] | list[bool], + a0: float, +) -> None: + """Mark points on periodic boundaries as fixed.""" + + if len(all_points) == 0: + return + + periodic_flags = [bool(value) for value in periodic] + periodic_mask = np.zeros(len(all_points), dtype=bool) + tolerance = max(0.05 * a0, BOUNDARY_FUZZ) + for axis, enabled in enumerate(periodic_flags): + if not enabled: + continue + periodic_mask |= np.isclose( + all_points[:, axis], geometry.bbox_min[axis], atol=tolerance, rtol=0.0 + ) + periodic_mask |= np.isclose( + all_points[:, axis], geometry.bbox_max[axis], atol=tolerance, rtol=0.0 ) - else: - generated_points = np.empty((0, geometry.dim), dtype=float) + states[periodic_mask] = STATE_FIXED - hint_points = [ - hint_points - for hint_points in geometry.piece_hints - if len(hint_points) > 0 - ] - if hint_points: - hint_block = _dedupe_points(np.vstack(hint_points)) - hint_block = filter_relevant(hint_block) - else: - hint_block = np.empty((0, geometry.dim), dtype=float) + +def _prepare_initial_points( + geometry: FemGeometry, + a0: float, + fixed_points: FloatArray, + mobile_points: FloatArray, + simply_points: FloatArray, + periodic: list[float] | list[bool], + rng: np.random.Generator, +) -> tuple[FloatArray, np.ndarray]: + """Prepare the initial point cloud and point-state array for relaxation.""" + + fixed_points, mobile_points, simply_points = _dedupe_fixed_mobile( + fixed_points, mobile_points, simply_points + ) + fixed_points = _filter_relevant_points(geometry, fixed_points) + mobile_points = _filter_relevant_points(geometry, mobile_points) + simply_points = _filter_relevant_points(geometry, simply_points) + generated_points = _select_generated_points( + geometry, + a0, + fixed_points, + mobile_points, + simply_points, + rng, + ) + hint_block = _collect_hint_points(geometry) all_points = np.vstack( ( @@ -199,27 +262,11 @@ def filter_relevant(points: FloatArray) -> FloatArray: ( np.full(len(fixed_points), STATE_FIXED, dtype=int), np.full(len(simply_points), STATE_SIMPLE, dtype=int), - np.full(len(mobile_points), STATE_MOBILE, dtype=int), + _classify_dynamic_states(geometry, mobile_points, a0), np.full(len(hint_block), STATE_FIXED, dtype=int), - np.full(len(generated_points), STATE_MOBILE, dtype=int), + _classify_dynamic_states(geometry, generated_points, a0), ) ) - if len(all_points) == 0: - return all_points, states - - periodic_flags = [bool(value) for value in periodic] - periodic_mask = np.zeros(len(all_points), dtype=bool) - tolerance = max(0.05 * a0, BOUNDARY_FUZZ) - for axis, enabled in enumerate(periodic_flags): - if not enabled: - continue - periodic_mask |= np.isclose( - all_points[:, axis], geometry.bbox_min[axis], atol=tolerance, rtol=0.0 - ) - periodic_mask |= np.isclose( - all_points[:, axis], geometry.bbox_max[axis], atol=tolerance, rtol=0.0 - ) - states[periodic_mask] = STATE_FIXED - + _apply_periodic_fixed_states(states, all_points, geometry, periodic, a0) return all_points, states diff --git a/src/nmesh/mesher/relaxation/topology.py b/src/nmesh/mesher/relaxation/topology.py index 7772a02..2b42674 100644 --- a/src/nmesh/mesher/relaxation/topology.py +++ b/src/nmesh/mesher/relaxation/topology.py @@ -11,7 +11,7 @@ from ...backend import RawMesh from ..periodic import build_periodic_groups -from ._constants import BOUNDARY_FUZZ, STATE_MOBILE +from ._constants import BOUNDARY_FUZZ, STATE_BOUNDARY, STATE_MOBILE from ._types import FloatArray from .geometry import FemGeometry @@ -94,6 +94,51 @@ def _region_volumes(region_ids: np.ndarray, measures: FloatArray) -> list[float] return [totals[region] for region in order] +def _regular_boundary_ratio(dim: int) -> float: + """Return the normalized volume-order ratio of an ideal regular simplex.""" + + if dim <= 0: + return 1.0 + return ((dim + 1) ** ((dim + 1) / 2.0)) / (math.factorial(dim) * (dim ** (dim / 2.0))) + + +def _simplex_volume_order_ratio(points: FloatArray, simplices: np.ndarray, dim: int) -> FloatArray: + """Return simplex volume divided by the local length scale raised to ``dim``.""" + + if len(simplices) == 0: + return np.empty(0, dtype=float) + + centroids = np.mean(points[simplices], axis=1) + offsets = points[simplices] - centroids[:, None, :] + max_radius = np.max(np.linalg.norm(offsets, axis=2), axis=1) + measures = _simplex_measures(points, simplices, dim) + order_scale = np.maximum(max_radius, BOUNDARY_FUZZ) ** dim + return measures / order_scale + + +def _classify_simplices_with_probes( + points: FloatArray, + simplices: np.ndarray, + geometry: FemGeometry, +) -> tuple[np.ndarray, np.ndarray]: + """Classify simplices using centroid and near-vertex probe points.""" + + if len(simplices) == 0: + return np.empty(0, dtype=int), np.empty(0, dtype=bool) + + centroids = np.mean(points[simplices], axis=1) + region_ids = geometry.classify_points(centroids) + offsets = points[simplices] - centroids[:, None, :] + probes = (centroids[:, None, :] + 0.9 * offsets).reshape(-1, geometry.dim) + probe_regions = geometry.classify_points(probes).reshape(len(simplices), geometry.dim + 1) + single_region = (not geometry.mesh_exterior) and (len(set(geometry.region_ids)) <= 1) + if single_region: + consistent = (region_ids >= 0) & np.all(probe_regions == region_ids[:, None], axis=1) + else: + consistent = (region_ids >= 0) & np.all(probe_regions >= 0, axis=1) + return region_ids, consistent + + def _triangulate_points(points: FloatArray, dim: int, states: np.ndarray | None = None) -> np.ndarray: """Triangulate the point cloud, retrying with light jitter for degenerate inputs.""" @@ -108,7 +153,11 @@ def _triangulate_points(points: FloatArray, dim: int, states: np.ndarray | None return Delaunay(points).simplices.astype(int, copy=False) except QhullError: jittered = np.array(points, copy=True) - movable = np.ones(len(points), dtype=bool) if states is None else states == STATE_MOBILE + movable = ( + np.ones(len(points), dtype=bool) + if states is None + else np.isin(states, [STATE_MOBILE, STATE_BOUNDARY]) + ) if np.any(movable): amplitudes = np.linspace(1.0e-9, 1.0e-8, np.count_nonzero(movable)) offsets = np.column_stack( @@ -128,18 +177,31 @@ def assemble_raw_mesh( points: FloatArray, geometry: FemGeometry, periodic: list[float] | list[bool], + *, + states: np.ndarray | None = None, + params: dict[str, Any] | None = None, ) -> RawMesh: """Assemble a ``RawMesh`` from relaxed points and geometry classification.""" coords = np.asarray(points, dtype=float) dim = geometry.dim - simplices = _triangulate_points(coords, dim) + simplices = _triangulate_points(coords, dim, states) if len(simplices) > 0: - centroids = np.mean(coords[simplices], axis=1) - region_ids = geometry.classify_points(centroids) + region_ids, probe_consistent = _classify_simplices_with_probes(coords, simplices, geometry) measures = _simplex_measures(coords, simplices, dim) - keep = (region_ids >= 0) & (measures > BOUNDARY_FUZZ) + boundary_mask = geometry.boundary_mask( + coords, + tolerance=max(BOUNDARY_FUZZ * 10.0, 1.0e-5), + ) + all_boundary = np.all(boundary_mask[simplices], axis=1) + boundary_ratio = _simplex_volume_order_ratio(coords, simplices, dim) + normalized_ratio = boundary_ratio / max(_regular_boundary_ratio(dim), BOUNDARY_FUZZ) + smallest_allowed_ratio = float( + (params or {}).get("controller_smallest_allowed_volume_ratio", 1.0) + ) + flat_boundary = all_boundary & (normalized_ratio < smallest_allowed_ratio) + keep = probe_consistent & (measures > BOUNDARY_FUZZ) & ~flat_boundary simplices = simplices[keep] region_ids = region_ids[keep] measures = measures[keep] diff --git a/tests/nmesh/test_integration.py b/tests/nmesh/test_integration.py new file mode 100644 index 0000000..48eb539 --- /dev/null +++ b/tests/nmesh/test_integration.py @@ -0,0 +1,58 @@ +import numpy as np + +import nmesh + + +def test_deterministic_meshing_uses_fixed_rng_seed(): + kwargs = dict( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[nmesh.Box([0.1, 0.1], [0.9, 0.9])], + a0=0.35, + ) + + mesh_a = nmesh.Mesh(**kwargs) + mesh_b = nmesh.Mesh(**kwargs) + + assert mesh_a.points == mesh_b.points + assert mesh_a.simplices == mesh_b.simplices + assert mesh_a.regions == mesh_b.regions + + +def test_periodic_scenario_completes_with_valid_topology(): + mesh = nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[nmesh.Box([0.25, 0.25], [0.75, 0.75])], + mesh_bounding_box=True, + periodic=[True, False], + a0=0.4, + ) + + assert len(mesh.points) > 0 + assert len(mesh.simplices) > 0 + assert len(mesh.periodic_point_indices) > 0 + assert all(len(group) >= 2 for group in mesh.periodic_point_indices) + assert all(region >= 0 for region in mesh.regions) + + +def test_hint_driven_scenario_completes_with_valid_topology(): + seed_mesh = nmesh.mesh_from_points_and_simplices( + points=[[0.2, 0.2], [0.8, 0.2], [0.5, 0.75]], + simplices_indices=[[0, 1, 2]], + simplices_regions=[5], + ) + + mesh = nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[], + mesh_bounding_box=True, + hints=[(seed_mesh, nmesh.Box([0.1, 0.1], [0.9, 0.9]))], + a0=0.4, + ) + + points = np.asarray(mesh.points, dtype=float) + assert len(mesh.points) > 0 + assert len(mesh.simplices) > 0 + assert all(mesh.points.count(point) == 1 for point in mesh.points) + assert np.any(np.all(np.isclose(points, [0.2, 0.2]), axis=1)) + assert np.any(np.all(np.isclose(points, [0.8, 0.2]), axis=1)) + assert np.any(np.all(np.isclose(points, [0.5, 0.75]), axis=1)) diff --git a/tests/nmesh/test_meshing_engine.py b/tests/nmesh/test_meshing_engine.py index 48ad9ee..1df9785 100644 --- a/tests/nmesh/test_meshing_engine.py +++ b/tests/nmesh/test_meshing_engine.py @@ -3,10 +3,173 @@ import numpy as np from nmesh import Box, Mesh -from nmesh.mesher.relaxation import _compile_density_function +from nmesh.mesher.relaxation import ( + _compile_density_function, + assemble_raw_mesh, + fem_geometry_from_bodies, +) +from nmesh.mesher.relaxation._constants import STATE_BOUNDARY, STATE_MOBILE +from nmesh.mesher.relaxation.forces import compute_forces +from nmesh.mesher.relaxation.seeding import _classify_dynamic_states class TestMeshingEngine(unittest.TestCase): + def test_boundary_state_detection_marks_surface_points(self): + geometry = fem_geometry_from_bodies( + (np.asarray([0.0, 0.0]), np.asarray([1.0, 1.0])), + [Box([0.0, 0.0], [1.0, 1.0]).obj], + [], + ) + + points = np.asarray( + [ + [0.0, 0.5], + [0.5, 0.5], + ], + dtype=float, + ) + + states = _classify_dynamic_states(geometry, points, a0=0.5) + + self.assertEqual(states.tolist(), [STATE_BOUNDARY, STATE_MOBILE]) + + def test_force_summary_suppresses_shape_force_at_boundary_corners(self): + geometry = fem_geometry_from_bodies( + (np.asarray([0.0, 0.0]), np.asarray([1.0, 1.0])), + [Box([0.0, 0.0], [1.0, 1.0]).obj], + [], + ) + + points = np.asarray( + [ + [0.0, 0.0], + [1.0, 0.0], + [1.0, 1.0], + [0.0, 1.0], + [0.6, 0.45], + ], + dtype=float, + ) + states = np.full(len(points), STATE_MOBILE, dtype=int) + summary = compute_forces( + points, + states, + geometry, + a0=0.5, + params={ + "controller_shape_force_scale": 0.3, + "controller_volume_force_scale": 0.0, + "controller_neigh_force_scale": 0.0, + "controller_irrel_elem_force_scale": 0.0, + "controller_sliver_correction": 1.0, + "controller_initial_settling_steps": 100, + }, + step=10, + ) + + self.assertGreaterEqual(summary.simplices.shape[0], 2) + self.assertTrue(np.allclose(summary.total_forces[:4], 0.0)) + self.assertGreater(np.linalg.norm(summary.total_forces[4]), 0.0) + self.assertGreater(summary.max_effective_force, 0.0) + + def test_boundary_effective_force_uses_tangential_projection(self): + geometry = fem_geometry_from_bodies( + (np.asarray([0.0, 0.0]), np.asarray([1.0, 1.0])), + [Box([0.0, 0.0], [1.0, 1.0]).obj], + [], + ) + + points = np.asarray( + [ + [0.0, 0.5], + [0.2, 0.4], + [0.2, 0.6], + ], + dtype=float, + ) + states = np.asarray([STATE_BOUNDARY, STATE_MOBILE, STATE_MOBILE], dtype=int) + summary = compute_forces( + points, + states, + geometry, + a0=0.4, + params={ + "controller_shape_force_scale": 0.0, + "controller_volume_force_scale": 0.0, + "controller_neigh_force_scale": 1.0, + "controller_irrel_elem_force_scale": 0.0, + "controller_sliver_correction": 1.0, + "controller_initial_settling_steps": 100, + }, + step=10, + ) + + self.assertGreater(abs(summary.total_forces[0][0]), 0.0) + self.assertAlmostEqual(summary.total_forces[0][1], 0.0, places=7) + self.assertAlmostEqual(summary.point_effective_force[0], 0.0, places=7) + + def test_geometry_projects_outside_point_back_to_boundary(self): + geometry = fem_geometry_from_bodies( + (np.asarray([0.0, 0.0]), np.asarray([1.0, 1.0])), + [Box([0.0, 0.0], [1.0, 1.0]).obj], + [], + ) + + projected = geometry.project_segment_to_domain( + np.asarray([0.2, 0.5], dtype=float), + np.asarray([-0.2, 0.5], dtype=float), + ) + + self.assertGreaterEqual(geometry.classify_points(projected.reshape(1, -1))[0], 0) + self.assertAlmostEqual(projected[0], 0.0, places=4) + + def test_assemble_raw_mesh_filters_flat_boundary_sliver(self): + geometry = fem_geometry_from_bodies( + (np.asarray([0.0, 0.0]), np.asarray([1.0, 0.01])), + [Box([0.0, 0.0], [1.0, 0.01]).obj], + [], + ) + + raw_mesh = assemble_raw_mesh( + np.asarray( + [ + [0.0, 0.0], + [1.0, 0.0], + [0.5, 0.01], + ], + dtype=float, + ), + geometry, + [False, False], + params={"controller_smallest_allowed_volume_ratio": 1.0}, + ) + + self.assertEqual(raw_mesh.simplices, []) + + def test_assemble_raw_mesh_keeps_regular_boundary_simplex(self): + height = np.sqrt(3.0) / 2.0 + geometry = fem_geometry_from_bodies( + (np.asarray([0.0, 0.0]), np.asarray([1.0, height])), + [Box([0.0, 0.0], [1.0, height]).obj], + [], + ) + + raw_mesh = assemble_raw_mesh( + np.asarray( + [ + [0.0, 0.0], + [1.0, 0.0], + [0.5, height], + ], + dtype=float, + ), + geometry, + [False, False], + params={"controller_smallest_allowed_volume_ratio": 1.0}, + ) + + self.assertEqual(len(raw_mesh.simplices), 1) + def test_density_string_translation_supports_c_style_blocks(self): density = _compile_density_function( """ From 5a80dc5383a12e137a003db0131a4d42702cea83 Mon Sep 17 00:00:00 2001 From: Emanuel Pituch <32016786+epituch@users.noreply.github.com> Date: Thu, 7 May 2026 19:58:54 -0700 Subject: [PATCH 3/3] WIP section 4 --- PARITY_COMPLETION.md | 53 +++ README.md | 12 +- src/nmesh/mesher/__init__.py | 1 + src/nmesh/mesher/meshing_parameters.py | 6 + src/nmesh/mesher/parity.py | 545 ++++++++++++++++++++++++ src/nmesh/mesher/periodic.py | 70 +-- src/nmesh/mesher/relaxation/density.py | 4 +- src/nmesh/mesher/relaxation/engine.py | 253 +++++++++-- src/nmesh/mesher/relaxation/finalize.py | 78 ++++ src/nmesh/mesher/relaxation/forces.py | 174 +++++--- src/nmesh/mesher/relaxation/geometry.py | 63 +++ src/nmesh/mesher/relaxation/jit.py | 6 +- src/nmesh/mesher/relaxation/recovery.py | 286 +++++++++++++ src/nmesh/mesher/relaxation/seeding.py | 210 +++++++-- src/nmesh/mesher/relaxation/topology.py | 42 +- tests/nmesh/test_integration.py | 64 +++ tests/nmesh/test_meshing_engine.py | 319 +++++++++++++- tests/nmesh/test_parity.py | 154 +++++++ tests/nmesh_test.py | 14 +- tools/nmesh_parity_compare.py | 247 +++++++++++ 20 files changed, 2396 insertions(+), 205 deletions(-) create mode 100644 PARITY_COMPLETION.md create mode 100644 src/nmesh/mesher/parity.py create mode 100644 src/nmesh/mesher/relaxation/finalize.py create mode 100644 src/nmesh/mesher/relaxation/recovery.py create mode 100644 tests/nmesh/test_parity.py create mode 100644 tools/nmesh_parity_compare.py diff --git a/PARITY_COMPLETION.md b/PARITY_COMPLETION.md new file mode 100644 index 0000000..ea928d8 --- /dev/null +++ b/PARITY_COMPLETION.md @@ -0,0 +1,53 @@ +# NMesh Section 4 Parity Completion + +**Status:** exact standalone behavioral parity required +**Branch:** `nmesh-section4` +**Reviewed:** 2026-05-04 +**Reference sources:** `nmag-src/src/mesh.ml`, `nmag-src/src/snippets.ml`, `nmag-src/interface/nmeshlib/lib1.py` + +## Position + +The Python 3 `nmesh` port must implement the full behavior of the original OCaml/Python2 mesher without requiring OCaml at runtime. Missing behavior is implementation work, not an accepted end state. + +The target is exact behavioral parity, not approximate parity. Topology decisions, region assignment, point-state transitions, add/delete scheduling, boundary recovery, periodic grouping, and file-visible mesh metadata must match the legacy implementation for the same inputs. Numeric tolerances are only allowed for floating-point coordinate comparisons, and each tolerance must be tied to a specific floating-point or triangulator tie-break reason. + +The port may use Python, NumPy, SciPy, Numba, Rust, C++, or C as needed. The public Python API should remain stable, and compiled helpers should be isolated behind that API. + +## Completed In This Section 4 Pass + +- Removed the Python-side 60-step cap; configured `controller_step_limit_max` is honored. +- Implemented legacy controller cadence: square-number add/delete checks, relaxed add/delete thresholds, topology-threshold retriangulation, minimum-step convergence, force-equilibrium stopping, and 50-step post-change settling. +- Aligned force metrics with legacy behavior: density-scaled movement, centroid density for simplex forces, tangent-projected boundary effective force, corner suppression, and OCaml-style Voronoi density correction. +- Classified force-time simplices with the same centroid, near-vertex probe, `Boundary` state, and raw volume-order ratio rules used by the legacy relevant/irrelevant topology split. +- Restricted irrelevant-element forces to mobile nodes, matching the active legacy path. +- Replaced deterministic grid seeding with the legacy random density estimator, D-lattice sphere-packing node estimate, rejection sampler, and RNG consumption order. +- Replaced midpoint point insertion with Gaussian insertion around the source point using local effective rod length. +- Ported active boundary recovery behavior from `mirror_simplices`: 2D mirrored points and 3D boundary-edge midpoint prevention points. +- Seeded paired fixed points on periodic outer-box faces, fixed only exact periodic boundary nodes, and canonicalized multi-axis periodic equivalence groups. +- Added final high-density moving-point cleanup, outside dynamic point cleanup, and final boundary snapping before final assembly. +- Normalized final simplex orientation to the legacy positive-volume convention. +- Fixed flat-boundary filtering so it follows the legacy `Boundary` state and raw volume-order ratio rules instead of deleting all geometrically boundary-adjacent multi-region simplices. +- Changed unsupported density snippets to fail loudly rather than silently falling back to density `1.0`. +- Added parity coverage for coarse adjacent pieces, concave difference domains, multi-axis periodic groups, final cleanup, and force-time relevant/irrelevant classification. +- Added canonical mesh signature helpers so parity fixtures compare topology and metadata exactly, with coordinate tolerance explicitly opt-in. +- Added a parity comparison runner (`tools/nmesh_parity_compare.py`) that generates modern scenario meshes and compares them with legacy `.nmesh` artifacts or a configured legacy runner command. +- Tightened runtime matching to avoid mesh-size-scaled boundary drift tolerance and rounded non-periodic periodic-group coordinates. + +## Validation Gate + +These are verification tasks, not accepted limitations: + +- Build reference fixtures from legacy examples and compare canonicalized topology, regions, surfaces, links, periodic groups, controller decisions, and mesh metadata exactly using `nmesh.mesher.parity`. +- Compare coordinates with the narrowest practical floating-point tolerance, documenting each tolerance and why exact binary equality is not the right assertion. +- Run the modernization scenario matrix through `tools/nmesh_parity_compare.py` and record metric-level acceptance thresholds for physically equivalent meshes when exact topology is not expected. +- Complete reference validation for periodic outer-box workflows on multiple periodic directions and complex periodic entities. +- Expand multi-piece/interface parity fixtures to include more interface-prevention cases. +- Profile large 3D meshes with legacy defaults and move any proven hot path to a compiled Rust, C++, or C helper behind the Python API. + +## Sign-Off Criteria + +- Legacy reference examples produce exact canonical mesh topology and metadata parity. +- Periodic, hint-driven, concave, and multi-piece examples match legacy behavior exactly, apart from explicitly documented floating-point coordinate tolerance. +- Add/delete, retriangulation, boundary recovery, and final cleanup decisions are covered by parity tests. +- Unsupported density syntax never changes mesh density silently. +- The full project test suite and parity fixtures pass in the project venv. diff --git a/README.md b/README.md index 914dee6..624e7f4 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,8 @@ This is an in-progress port of [nmag](https://github.com/nmag-project/nmag-src) that updates old dependencies and packages nmag as a standalone python 3 module. +The Python 3 port is intended to be standalone and must not require OCaml at runtime. See [PARITY_COMPLETION.md](PARITY_COMPLETION.md) for the Section 4 parity record and validation gate. + ### Installation * Create a virtual environment: `python -m venv venv` * Activate the virtual environment: @@ -12,13 +14,7 @@ This is an in-progress port of [nmag](https://github.com/nmag-project/nmag-src) ### Testing * To install the optional test dependencies run `pip install -e ".[test]"` * Once the project has been initialized simply run `pytest` to run your tests +* To generate modern nmesh parity scenario outputs run `python tools/nmesh_parity_compare.py --output-dir parity/new` +* To compare against legacy `.nmesh` artifacts run `python tools/nmesh_parity_compare.py --legacy-dir parity/legacy` * To view coverage open `htmlcov/index.html` in your browser * [Specifying what tests to run](https://docs.pytest.org/en/latest/how-to/usage.html#specifying-which-tests-to-run) - -### OCaml -* Install latest Ocaml/Opam version 4 https://ocaml.org/install - * Recommended to run `opam init` and setup a switch with version 4 -* Install ocaml-in-python `opam install ocaml-in-python` -* Register the package in python ``pip install --editable "`opam var ocaml-in-python:lib`"`` -* Tell python where to look for the OCaml library `export OCAMLPATH=${DUNE_DIR}/_build/install/default/lib` where DUNE_DIR is the directory of the dune project -* Run this to explictly activate the Opam switch `eval $(opam env)` diff --git a/src/nmesh/mesher/__init__.py b/src/nmesh/mesher/__init__.py index 579ddf5..ed2e173 100644 --- a/src/nmesh/mesher/__init__.py +++ b/src/nmesh/mesher/__init__.py @@ -1,3 +1,4 @@ from .meshing_parameters import * from .driver import * +from .parity import * from .relaxation import * diff --git a/src/nmesh/mesher/meshing_parameters.py b/src/nmesh/mesher/meshing_parameters.py index d366013..664c26c 100644 --- a/src/nmesh/mesher/meshing_parameters.py +++ b/src/nmesh/mesher/meshing_parameters.py @@ -341,6 +341,9 @@ def to_mesher_config(self, dim): continue resolved[spec.internal_name] = spec.cast(value) + for key, value in self._params.items(): + resolved.setdefault(key, value) + return resolved def apply_to_mesher(self, mesher, dim): @@ -354,6 +357,9 @@ def apply_to_mesher(self, mesher, dim): continue mesher["parameters"][spec.internal_name] = spec.cast(value) + for key, value in self._params.items(): + mesher["parameters"].setdefault(key, value) + return mesher def _set_parameter(self, name, value): diff --git a/src/nmesh/mesher/parity.py b/src/nmesh/mesher/parity.py new file mode 100644 index 0000000..de15198 --- /dev/null +++ b/src/nmesh/mesher/parity.py @@ -0,0 +1,545 @@ +"""Canonical mesh signatures for exact parity comparisons.""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path +from typing import TypeAlias + +import numpy as np + +from ..backend import RawMesh + +CoordinateKey: TypeAlias = tuple[str | int, ...] + +__all__ = [ + "CanonicalMeshSignature", + "MeshMetricComparison", + "MeshMetricSummary", + "assert_canonical_mesh_equal", + "canonical_mesh_signature", + "compare_mesh_metrics", + "mesh_metric_summary", + "read_ascii_nmesh", +] + + +@dataclass(frozen=True, slots=True) +class CanonicalMeshSignature: + """Order-independent mesh signature for legacy parity checks. + + Topology and metadata are compared exactly after canonical point reindexing. + Coordinate tolerance is opt-in and only affects coordinate keys; leaving it + as ``None`` uses exact binary floating-point identity. + """ + + dim: int + point_keys: tuple[CoordinateKey, ...] + point_regions: tuple[tuple[int, ...], ...] + simplices: tuple[tuple[int, tuple[int, ...]], ...] + surfaces: tuple[tuple[int, ...], ...] + links: tuple[tuple[int, int], ...] + periodic_groups: tuple[tuple[int, ...], ...] + region_volumes: tuple[float, ...] + + +@dataclass(frozen=True, slots=True) +class MeshMetricSummary: + """Geometry-oriented mesh summary for modernization parity checks.""" + + dim: int + point_count: int + simplex_count: int + surface_count: int + link_count: int + periodic_group_count: int + region_counts: tuple[tuple[int, int], ...] + region_volumes: tuple[tuple[int, float], ...] + bbox_min: tuple[float, ...] + bbox_max: tuple[float, ...] + simplex_measure_min: float + simplex_measure_mean: float + simplex_measure_max: float + edge_length_min: float + edge_length_mean: float + edge_length_max: float + + +@dataclass(frozen=True, slots=True) +class MeshMetricComparison: + """Result of a metric-level comparison between two meshes.""" + + actual: MeshMetricSummary + expected: MeshMetricSummary + failures: tuple[str, ...] + + @property + def passed(self) -> bool: + """Return whether all configured metric checks passed.""" + + return not self.failures + + +def canonical_mesh_signature( + raw_mesh: RawMesh, + *, + coordinate_tolerance: float | None = None, +) -> CanonicalMeshSignature: + """Return a canonical signature suitable for exact parity assertions.""" + + point_keys = [_coordinate_key(point, coordinate_tolerance) for point in raw_mesh.points] + canonical_order = sorted(range(len(point_keys)), key=point_keys.__getitem__) + canonical_index = {old_index: new_index for new_index, old_index in enumerate(canonical_order)} + _require_unique_points(point_keys) + + ordered_point_keys = tuple(point_keys[index] for index in canonical_order) + ordered_point_regions = tuple( + tuple(sorted(map(int, _point_regions_at(raw_mesh, index)))) + for index in canonical_order + ) + canonical_simplices = tuple( + sorted( + ( + int(region), + tuple(canonical_index[int(point_index)] for point_index in simplex), + ) + for region, simplex in zip(raw_mesh.regions, raw_mesh.simplices) + ) + ) + canonical_surfaces = tuple( + sorted( + tuple(sorted(canonical_index[int(point_index)] for point_index in surface)) + for surface in raw_mesh.surfaces + ) + ) + canonical_links = tuple( + sorted( + tuple(sorted((canonical_index[int(left)], canonical_index[int(right)]))) + for left, right in raw_mesh.links + ) + ) + canonical_periodic = tuple( + sorted( + tuple(sorted(canonical_index[int(point_index)] for point_index in group)) + for group in raw_mesh.periodic_point_indices + ) + ) + return CanonicalMeshSignature( + dim=int(raw_mesh.dim), + point_keys=ordered_point_keys, + point_regions=ordered_point_regions, + simplices=canonical_simplices, + surfaces=canonical_surfaces, + links=canonical_links, + periodic_groups=canonical_periodic, + region_volumes=tuple(float(volume) for volume in raw_mesh.region_volumes), + ) + + +def assert_canonical_mesh_equal( + actual: RawMesh, + expected: RawMesh, + *, + coordinate_tolerance: float | None = None, +) -> None: + """Assert exact canonical parity between two meshes.""" + + actual_signature = canonical_mesh_signature( + actual, + coordinate_tolerance=coordinate_tolerance, + ) + expected_signature = canonical_mesh_signature( + expected, + coordinate_tolerance=coordinate_tolerance, + ) + if actual_signature != expected_signature: + raise AssertionError( + "Canonical mesh signatures differ:\n" + f"actual={actual_signature!r}\n" + f"expected={expected_signature!r}" + ) + + +def _coordinate_key( + point: list[float], + coordinate_tolerance: float | None, +) -> CoordinateKey: + """Return an exact or explicitly quantized coordinate key.""" + + if coordinate_tolerance is None: + return tuple(float(value).hex() for value in point) + if coordinate_tolerance <= 0.0: + raise ValueError("coordinate_tolerance must be positive when provided") + return tuple(int(round(float(value) / coordinate_tolerance)) for value in point) + + +def _require_unique_points(point_keys: list[CoordinateKey]) -> None: + """Reject signatures where canonical point identity is ambiguous.""" + + if len(set(point_keys)) != len(point_keys): + raise ValueError( + "Canonical mesh signature requires unique point coordinates at the " + "selected coordinate tolerance" + ) + + +def _point_regions_at(raw_mesh: RawMesh, index: int) -> list[int]: + """Return point-region membership for a point, tolerating absent metadata.""" + + if index < len(raw_mesh.point_regions): + return raw_mesh.point_regions[index] + return [] + + +def read_ascii_nmesh(path: str | Path) -> RawMesh: + """Read the legacy ASCII ``.nmesh``/PYFEM mesh format into ``RawMesh``.""" + + source = Path(path) + lines = [ + line.strip() + for line in source.read_text(encoding="utf-8").splitlines() + if line.strip() and not line.lstrip().startswith("#") + ] + cursor = 0 + point_count = int(lines[cursor]) + cursor += 1 + points = [_parse_float_row(lines[cursor + index]) for index in range(point_count)] + cursor += point_count + dim = len(points[0]) if points else 0 + + simplex_count = int(lines[cursor]) + cursor += 1 + simplices: list[list[int]] = [] + regions: list[int] = [] + for _ in range(simplex_count): + values = _parse_int_row(lines[cursor]) + cursor += 1 + regions.append(values[0]) + simplices.append(values[1:]) + + surface_count = int(lines[cursor]) + cursor += 1 + surfaces: list[list[int]] = [] + for _ in range(surface_count): + values = _parse_int_row(lines[cursor]) + cursor += 1 + surfaces.append(values[-dim:] if dim > 0 else []) + + periodic_groups: list[list[int]] = [] + if cursor < len(lines): + periodic_count = int(lines[cursor]) + cursor += 1 + for _ in range(periodic_count): + values = _parse_int_row(lines[cursor]) + cursor += 1 + periodic_groups.append(values[1:]) + + return RawMesh( + points=points, + simplices=simplices, + regions=regions, + point_regions=_build_point_regions(len(points), simplices, regions), + surfaces=surfaces, + links=_build_links(simplices), + region_volumes=_region_volumes(points, simplices, regions, dim), + periodic_point_indices=periodic_groups, + permutation=list(range(len(points))), + dim=dim, + ) + + +def mesh_metric_summary(raw_mesh: RawMesh) -> MeshMetricSummary: + """Return numeric mesh metrics useful for legacy-vs-modern comparisons.""" + + points = np.asarray(raw_mesh.points, dtype=float) + dim = int(raw_mesh.dim) + if points.size == 0: + points = np.empty((0, dim), dtype=float) + elif points.ndim == 1: + points = points.reshape(1, dim) + + simplices = np.asarray(raw_mesh.simplices, dtype=int) + simplex_measures = _simplex_measures(points, simplices, dim) + edge_lengths = _edge_lengths(points, raw_mesh.links) + region_counts = tuple( + sorted( + (int(region), int(sum(1 for value in raw_mesh.regions if int(value) == int(region)))) + for region in set(raw_mesh.regions) + ) + ) + region_volumes = tuple( + sorted(zip(sorted(set(map(int, raw_mesh.regions))), map(float, raw_mesh.region_volumes))) + ) + return MeshMetricSummary( + dim=dim, + point_count=len(raw_mesh.points), + simplex_count=len(raw_mesh.simplices), + surface_count=len(raw_mesh.surfaces), + link_count=len(raw_mesh.links), + periodic_group_count=len(raw_mesh.periodic_point_indices), + region_counts=region_counts, + region_volumes=region_volumes, + bbox_min=tuple(np.min(points, axis=0).tolist()) if len(points) else (), + bbox_max=tuple(np.max(points, axis=0).tolist()) if len(points) else (), + simplex_measure_min=_safe_min(simplex_measures), + simplex_measure_mean=_safe_mean(simplex_measures), + simplex_measure_max=_safe_max(simplex_measures), + edge_length_min=_safe_min(edge_lengths), + edge_length_mean=_safe_mean(edge_lengths), + edge_length_max=_safe_max(edge_lengths), + ) + + +def compare_mesh_metrics( + actual: RawMesh, + expected: RawMesh, + *, + count_relative_tolerance: float = 0.25, + volume_relative_tolerance: float = 0.05, + length_relative_tolerance: float = 0.10, + absolute_tolerance: float = 1.0e-9, +) -> MeshMetricComparison: + """Compare two meshes with documented metric-level tolerances.""" + + actual_summary = mesh_metric_summary(actual) + expected_summary = mesh_metric_summary(expected) + failures: list[str] = [] + _compare_exact("dim", actual_summary.dim, expected_summary.dim, failures) + _compare_relative_count( + "point_count", + actual_summary.point_count, + expected_summary.point_count, + count_relative_tolerance, + failures, + ) + _compare_relative_count( + "simplex_count", + actual_summary.simplex_count, + expected_summary.simplex_count, + count_relative_tolerance, + failures, + ) + _compare_relative_count( + "surface_count", + actual_summary.surface_count, + expected_summary.surface_count, + count_relative_tolerance, + failures, + ) + _compare_exact("region ids", _region_ids(actual_summary), _region_ids(expected_summary), failures) + _compare_exact( + "periodic_group_count", + actual_summary.periodic_group_count, + expected_summary.periodic_group_count, + failures, + ) + _compare_vector( + "bbox_min", + actual_summary.bbox_min, + expected_summary.bbox_min, + length_relative_tolerance, + absolute_tolerance, + failures, + ) + _compare_vector( + "bbox_max", + actual_summary.bbox_max, + expected_summary.bbox_max, + length_relative_tolerance, + absolute_tolerance, + failures, + ) + _compare_region_volumes( + actual_summary.region_volumes, + expected_summary.region_volumes, + volume_relative_tolerance, + absolute_tolerance, + failures, + ) + for field_name in ( + "simplex_measure_min", + "simplex_measure_mean", + "simplex_measure_max", + "edge_length_min", + "edge_length_mean", + "edge_length_max", + ): + _compare_float( + field_name, + float(getattr(actual_summary, field_name)), + float(getattr(expected_summary, field_name)), + length_relative_tolerance, + absolute_tolerance, + failures, + ) + return MeshMetricComparison( + actual=actual_summary, + expected=expected_summary, + failures=tuple(failures), + ) + + +def _parse_float_row(line: str) -> list[float]: + return [float(value) for value in line.split()] + + +def _parse_int_row(line: str) -> list[int]: + return [int(value) for value in line.split()] + + +def _build_point_regions( + point_count: int, + simplices: list[list[int]], + regions: list[int], +) -> list[list[int]]: + memberships: list[set[int]] = [set() for _ in range(point_count)] + for simplex, region in zip(simplices, regions): + for point_index in simplex: + memberships[int(point_index)].add(int(region)) + return [sorted(group) for group in memberships] + + +def _build_links(simplices: list[list[int]]) -> list[tuple[int, int]]: + links: set[tuple[int, int]] = set() + for simplex in simplices: + for left_index, left in enumerate(simplex): + for right in simplex[left_index + 1 :]: + a = int(left) + b = int(right) + links.add((a, b) if a <= b else (b, a)) + return sorted(links) + + +def _region_volumes( + points: list[list[float]], + simplices: list[list[int]], + regions: list[int], + dim: int, +) -> list[float]: + if not simplices: + return [] + coords = np.asarray(points, dtype=float) + simplex_array = np.asarray(simplices, dtype=int) + measures = _simplex_measures(coords, simplex_array, dim) + totals = {int(region): 0.0 for region in regions} + for region, measure in zip(regions, measures): + totals[int(region)] += float(measure) + return [totals[region] for region in sorted(totals)] + + +def _simplex_measures(points: np.ndarray, simplices: np.ndarray, dim: int) -> np.ndarray: + if len(simplices) == 0: + return np.empty(0, dtype=float) + if dim == 1: + return np.abs(points[simplices[:, 1], 0] - points[simplices[:, 0], 0]) + matrices = points[simplices[:, 1:]] - points[simplices[:, [0]]] + return np.abs(np.linalg.det(matrices)) / float(_factorial(dim)) + + +def _edge_lengths(points: np.ndarray, links: list[tuple[int, int]]) -> np.ndarray: + if not links: + return np.empty(0, dtype=float) + return np.asarray( + [np.linalg.norm(points[int(left)] - points[int(right)]) for left, right in links], + dtype=float, + ) + + +def _factorial(value: int) -> int: + result = 1 + for item in range(2, value + 1): + result *= item + return result + + +def _safe_min(values: np.ndarray) -> float: + return float(np.min(values)) if len(values) else 0.0 + + +def _safe_mean(values: np.ndarray) -> float: + return float(np.mean(values)) if len(values) else 0.0 + + +def _safe_max(values: np.ndarray) -> float: + return float(np.max(values)) if len(values) else 0.0 + + +def _compare_exact(name: str, actual: object, expected: object, failures: list[str]) -> None: + if actual != expected: + failures.append(f"{name}: actual={actual!r}, expected={expected!r}") + + +def _compare_relative_count( + name: str, + actual: int, + expected: int, + relative_tolerance: float, + failures: list[str], +) -> None: + allowed = max(1.0, abs(float(expected)) * relative_tolerance) + if abs(float(actual - expected)) > allowed: + failures.append( + f"{name}: actual={actual}, expected={expected}, allowed_delta={allowed:g}" + ) + + +def _compare_vector( + name: str, + actual: tuple[float, ...], + expected: tuple[float, ...], + relative_tolerance: float, + absolute_tolerance: float, + failures: list[str], +) -> None: + if len(actual) != len(expected): + failures.append(f"{name}: actual={actual!r}, expected={expected!r}") + return + for index, (actual_value, expected_value) in enumerate(zip(actual, expected)): + _compare_float( + f"{name}[{index}]", + actual_value, + expected_value, + relative_tolerance, + absolute_tolerance, + failures, + ) + + +def _compare_region_volumes( + actual: tuple[tuple[int, float], ...], + expected: tuple[tuple[int, float], ...], + relative_tolerance: float, + absolute_tolerance: float, + failures: list[str], +) -> None: + if tuple(region for region, _volume in actual) != tuple(region for region, _volume in expected): + failures.append(f"region volumes ids: actual={actual!r}, expected={expected!r}") + return + for (region, actual_volume), (_expected_region, expected_volume) in zip(actual, expected): + _compare_float( + f"region_volumes[{region}]", + actual_volume, + expected_volume, + relative_tolerance, + absolute_tolerance, + failures, + ) + + +def _compare_float( + name: str, + actual: float, + expected: float, + relative_tolerance: float, + absolute_tolerance: float, + failures: list[str], +) -> None: + allowed = max(absolute_tolerance, abs(expected) * relative_tolerance) + if abs(actual - expected) > allowed: + failures.append( + f"{name}: actual={actual:g}, expected={expected:g}, allowed_delta={allowed:g}" + ) + + +def _region_ids(summary: MeshMetricSummary) -> tuple[int, ...]: + return tuple(region for region, _count in summary.region_counts) diff --git a/src/nmesh/mesher/periodic.py b/src/nmesh/mesher/periodic.py index a4ecc04..0d89739 100644 --- a/src/nmesh/mesher/periodic.py +++ b/src/nmesh/mesher/periodic.py @@ -1,19 +1,13 @@ -from __future__ import annotations - """Helpers for building periodic-equivalence groups from mesh coordinates.""" +from __future__ import annotations + from collections import defaultdict from collections.abc import Iterable import numpy as np -def _point_key(values: np.ndarray, decimals: int = 8) -> tuple[float, ...]: - """Return a rounded tuple key suitable for coordinate-based lookups.""" - - return tuple(np.round(values.astype(float, copy=False), decimals=decimals).tolist()) - - def merge_periodic_groups(groups: Iterable[Iterable[int]]) -> list[list[int]]: """Merge overlapping periodic point groups into disjoint sorted groups.""" @@ -54,28 +48,50 @@ def build_periodic_groups( return [] periodic_flags = [bool(value) for value in periodic_axes] - raw_groups: list[list[int]] = [] - - for axis, enabled in enumerate(periodic_flags): - if not enabled: - continue + if not any(periodic_flags): + return [] - min_mask = np.isclose(points[:, axis], bbox_min[axis], atol=tolerance, rtol=0.0) - max_mask = np.isclose(points[:, axis], bbox_max[axis], atol=tolerance, rtol=0.0) - if not np.any(min_mask) or not np.any(max_mask): + grouped: defaultdict[tuple[tuple[str, float | int], ...], list[int]] = defaultdict(list) + for index, point in enumerate(points): + key = _periodic_equivalence_key( + point, + bbox_min, + bbox_max, + periodic_flags, + tolerance, + ) + if key is None: continue + grouped[key].append(int(index)) - other_axes = [other for other in range(points.shape[1]) if other != axis] - min_points = np.flatnonzero(min_mask) - max_points = np.flatnonzero(max_mask) + return [sorted(group) for group in grouped.values() if len(group) >= 2] - max_lookup: defaultdict[tuple[float, ...], list[int]] = defaultdict(list) - for index in max_points: - max_lookup[_point_key(points[index, other_axes])].append(int(index)) - for index in min_points: - key = _point_key(points[index, other_axes]) - for partner in max_lookup.get(key, []): - raw_groups.append([int(index), int(partner)]) +def _periodic_equivalence_key( + point: np.ndarray, + bbox_min: np.ndarray, + bbox_max: np.ndarray, + periodic_flags: list[bool], + tolerance: float, +) -> tuple[tuple[str, float | int], ...] | None: + """Return a canonical key for all periodic copies of one boundary point.""" + + key: list[tuple[str, float | int]] = [] + touches_periodic_boundary = False + for axis, value in enumerate(point): + coord = float(value) + if not periodic_flags[axis]: + key.append(("coord", coord)) + continue - return merge_periodic_groups(raw_groups) + on_min = np.isclose(coord, bbox_min[axis], atol=tolerance, rtol=0.0) + on_max = np.isclose(coord, bbox_max[axis], atol=tolerance, rtol=0.0) + if on_min or on_max: + key.append(("periodic", axis)) + touches_periodic_boundary = True + else: + key.append(("coord", coord)) + + if not touches_periodic_boundary: + return None + return tuple(key) diff --git a/src/nmesh/mesher/relaxation/density.py b/src/nmesh/mesher/relaxation/density.py index 6f9d4bf..86ca2c1 100644 --- a/src/nmesh/mesher/relaxation/density.py +++ b/src/nmesh/mesher/relaxation/density.py @@ -162,5 +162,5 @@ def _compile_density_function(density: str | DensityFunction | None) -> DensityF if expression_density is not None: return expression_density return _compile_script_density(source) - except Exception: - return lambda _point: 1.0 + except Exception as exc: + raise ValueError(f"Could not compile density function: {source!r}") from exc diff --git a/src/nmesh/mesher/relaxation/engine.py b/src/nmesh/mesher/relaxation/engine.py index e468f58..7bc8e58 100644 --- a/src/nmesh/mesher/relaxation/engine.py +++ b/src/nmesh/mesher/relaxation/engine.py @@ -14,8 +14,10 @@ from ..meshing_parameters import PointFate, default_handle_point_density_fun from ._constants import BOUNDARY_FUZZ, DEFAULT_RNG_SEED, STATE_BOUNDARY, STATE_FIXED, STATE_MOBILE, STATE_SIMPLE from ._types import DensityFunction, EngineResult, FloatArray +from .finalize import snap_final_boundary_points from .forces import ForceSummary, compute_forces from .geometry import FemGeometry, fem_geometry_from_bodies +from .recovery import mirror_surface_recovery_points from .seeding import _as_float_array, _classify_dynamic_states, _dedupe_points, _prepare_initial_points from .topology import _callback_mesh_info, assemble_raw_mesh @@ -52,11 +54,17 @@ def __init__( simply_points, self.periodic, self.rng, + self.params, ) self.step = 0 self.finished = False + self.last_addition_deletion_step = 0 self.last_max_displacement = math.inf self.last_max_effective_force = math.inf + self.last_point_density = np.zeros(len(self.points), dtype=float) + self.max_rel_movement_since_last_triangulation = 0.0 + self.current_simplices: np.ndarray | None = None + self.points_at_last_triangulation = np.array(self.points, copy=True) self._refresh_boundary_states() self.cached_raw_mesh = assemble_raw_mesh( self.points, @@ -68,9 +76,9 @@ def __init__( @property def max_steps(self) -> int: - """Return the capped step budget for this Python port.""" + """Return the configured maximum relaxation-step budget.""" - return min(int(self.params.get("controller_step_limit_max", 1000)), 60) + return int(self.params.get("controller_step_limit_max", 1000)) @property def tolerated_rel_move(self) -> float: @@ -94,32 +102,67 @@ def max_time_step(self) -> float: def boundary_drift_tolerance(self) -> float: """Return the tolerated increase in boundary distance for boundary points. - We allow a small positive tolerance here rather than requiring a strict - monotonic decrease on every step. That matches the practical numerical - behavior of the Python port better and avoids jitter from floating-point - noise when a point is already very close to the surface. + Legacy boundary points are only allowed to move when they do not drift + farther from the boundary. The Python port keeps only the global + floating-point boundary fuzz here; it does not scale this tolerance with + ``a0`` because that would turn a numerical guard into different meshing + behavior. """ - return max(BOUNDARY_FUZZ, 0.02 * self.a0) + return BOUNDARY_FUZZ @property def min_equilibrium_steps(self) -> int: """Return the minimum number of steps before equilibrium can stop the loop.""" - configured = int(self.params.get("controller_step_limit_min", 500)) - return min(configured, min(4, self.max_steps)) + return int(self.params.get("controller_step_limit_min", 500)) + + @property + def post_change_settling_steps(self) -> int: + """Return the legacy settling window after point add/delete checks.""" + + return 50 - def _rebuild_mesh(self) -> None: + def _rebuild_mesh(self, *, final: bool = False) -> None: """Refresh the cached ``RawMesh`` snapshot from the current point cloud.""" + points, states = self._final_output_points_and_states() if final else (self.points, self.states) self.cached_raw_mesh = assemble_raw_mesh( - self.points, + points, self.geometry, self.periodic, - states=self.states, + states=states, params=self.params, ) + def _final_output_points_and_states(self) -> tuple[FloatArray, np.ndarray]: + """Return the final point cloud after the legacy high-density cleanup.""" + + if len(self.last_point_density) != len(self.points): + return self.points, self.states + dynamic_mask = np.isin(self.states, [STATE_MOBILE, STATE_BOUNDARY]) + keep = (~dynamic_mask) | (self.last_point_density < 100.0) + if np.all(keep): + points = self.points + states = self.states + else: + points = self.points[keep] + states = self.states[keep] + dynamic_mask = np.isin(states, [STATE_MOBILE, STATE_BOUNDARY]) + domain_mask = self.geometry.classify_points(points) >= 0 + keep_domain = (~dynamic_mask) | domain_mask + if not np.all(keep_domain): + points = points[keep_domain] + states = states[keep_domain] + return snap_final_boundary_points(points, states, self.geometry, self.a0, self.params) + + def _mark_topology_stale(self) -> None: + """Request a fresh Delaunay triangulation before the next force step.""" + + self.current_simplices = None + self.max_rel_movement_since_last_triangulation = 0.0 + self.points_at_last_triangulation = np.array(self.points, copy=True) + def _refresh_boundary_states(self) -> None: """Refresh dynamic points between mobile and boundary states.""" @@ -132,19 +175,110 @@ def _refresh_boundary_states(self) -> None: def _effective_time_step(self, force_summary: ForceSummary) -> float: """Compute the controller-like time step for the current force field.""" - settling_steps = int(self.params.get("controller_initial_settling_steps", 100)) - weight_fun = self.params.get("initial_relaxation_weight_fun") - if callable(weight_fun): - relaxation_weight = float(weight_fun(self.step, settling_steps, 0.0, 1.0)) - else: - relaxation_weight = min(1.0, float(self.step) / max(float(settling_steps), 1.0)) - capped_max_time_step = relaxation_weight * self.max_time_step + freedom = float(self.params.get("controller_movement_max_freedom", 3.0)) + movement_weight = self._initial_relaxation_weight(freedom, 1.0) + capped_max_time_step = movement_weight * self.max_time_step if force_summary.max_effective_force <= BOUNDARY_FUZZ: return capped_max_time_step return min( capped_max_time_step, - (relaxation_weight * self.time_step_scale) / force_summary.max_effective_force, + (movement_weight * self.time_step_scale) / force_summary.max_effective_force, + ) + + def _initial_relaxation_weight(self, init_val: float, final_val: float) -> float: + """Evaluate the legacy initial-relaxation ramp for this step.""" + + settling_steps = int(self.params.get("controller_initial_settling_steps", 100)) + weight_fun = self.params.get("initial_relaxation_weight_fun") + if callable(weight_fun): + return float(weight_fun(self.step, settling_steps, init_val, final_val)) + fraction = min(1.0, float(self.step) / max(float(settling_steps), 1.0)) + return init_val + (final_val - init_val) * fraction + + def _density_thresholds(self) -> tuple[float, float]: + """Return add/delete thresholds with the legacy initial relaxation ramp.""" + + freedom = float(self.params.get("controller_movement_max_freedom", 3.0)) + thresh_add = ( + self._initial_relaxation_weight(-0.1 * freedom, 0.0) + + float(self.params.get("controller_thresh_add", 1.0)) + ) + thresh_del = ( + self._initial_relaxation_weight(0.1 * freedom, 0.0) + + float(self.params.get("controller_thresh_del", 2.0)) ) + return thresh_add, thresh_del + + @staticmethod + def _is_positive_square_number(value: int) -> bool: + """Return whether ``value`` is a positive perfect square.""" + + if value <= 0: + return False + root = int(math.sqrt(value) + 0.5) + return root * root == value + + def _should_attempt_point_change(self) -> bool: + """Return whether the legacy controller schedules point fate changes.""" + + return ( + self.step < self.max_steps + and self._is_positive_square_number(self.step - 10) + ) + + def _step_limit_reached(self) -> bool: + """Return whether the legacy max-step and post-change settling rules are met.""" + + return ( + self.step >= self.max_steps + and self.step >= self.last_addition_deletion_step + self.post_change_settling_steps + ) + + def _topology_threshold(self) -> float: + """Return the relaxed legacy threshold for topology refresh.""" + + freedom = float(self.params.get("controller_movement_max_freedom", 3.0)) + threshold = float(self.params.get("controller_topology_threshold", 0.2)) + return self._initial_relaxation_weight(freedom, 1.0) * threshold + + def _record_triangulation(self, force_summary: ForceSummary) -> None: + """Cache the current topology after a force calculation builds it.""" + + if self.current_simplices is None: + self.current_simplices = np.asarray(force_summary.simplices, dtype=int) + self.points_at_last_triangulation = np.array(self.points, copy=True) + self.max_rel_movement_since_last_triangulation = 0.0 + + def _update_topology_movement(self) -> None: + """Track movement since the last triangulation in OCaml controller units.""" + + if len(self.points) != len(self.points_at_last_triangulation): + self._mark_topology_stale() + return + + if len(self.points) == 0: + self.max_rel_movement_since_last_triangulation = 0.0 + return + + displacements = np.linalg.norm(self.points - self.points_at_last_triangulation, axis=1) + density_scale = np.asarray( + [ + self.geometry.density_at(point) ** (1.0 / max(self.geometry.dim, 1)) + for point in self.points + ], + dtype=float, + ) + scaled = displacements * density_scale / max(self.a0, BOUNDARY_FUZZ) + self.max_rel_movement_since_last_triangulation = max( + self.max_rel_movement_since_last_triangulation, + float(np.max(scaled, initial=0.0)), + ) + + def _refresh_topology_if_needed(self) -> None: + """Invalidate cached topology when movement exceeds the legacy threshold.""" + + if self.max_rel_movement_since_last_triangulation > self._topology_threshold(): + self._mark_topology_stale() def _attempt_add_delete_points(self, force_summary: ForceSummary) -> None: """Apply the mesher density heuristic to add or remove mobile points.""" @@ -153,6 +287,15 @@ def _attempt_add_delete_points(self, force_summary: ForceSummary) -> None: return additions, removals = self._evaluate_point_densities(force_summary) + recovery_points = mirror_surface_recovery_points( + self.points, + self.states, + self.geometry, + self.a0, + force_summary.simplices, + ) + if len(recovery_points) > 0: + additions.extend(recovery_points) self._apply_point_changes(additions, removals) self._refresh_boundary_states() @@ -163,8 +306,7 @@ def _evaluate_point_densities( """Decide which points to add or remove from the current cloud.""" handler = self.params.get("handle_point_density_fun", default_handle_point_density_fun) - thresh_add = float(self.params.get("controller_thresh_add", 1.0)) - thresh_del = float(self.params.get("controller_thresh_del", 2.0)) + thresh_add, thresh_del = self._density_thresholds() additions: list[FloatArray] = [] removals: list[int] = [] @@ -177,22 +319,30 @@ def _evaluate_point_densities( continue point = self.points[index] - neigh_coords = self.points[np.asarray(neigh, dtype=int)] - distances = np.linalg.norm(neigh_coords - point, axis=1) avg_density = float(force_summary.point_density[index]) avg_force = float(force_summary.point_average_force[index]) fate = handler(self.rng, (avg_density, avg_force), thresh_add, thresh_del) if fate == PointFate.ADD_ANOTHER: - farthest = neigh_coords[int(np.argmax(distances))] - candidate = 0.5 * (point + farthest) - if self.geometry.classify_points(candidate.reshape(1, -1))[0] >= 0: - additions.append(candidate) + additions.append(self._random_point_close_to(point)) elif fate == PointFate.DELETE and len(self.points) - len(removals) > self.geometry.dim + 1: removals.append(index) return additions, removals + def _random_point_close_to(self, point: FloatArray) -> FloatArray: + """Return a Gaussian insertion candidate using the legacy rod-length scale.""" + + density_here = self.geometry.density_at(point) + effective_rod_length = self.a0 * (density_here ** (-1.0 / max(self.geometry.dim, 1))) + candidate = np.asarray( + self.rng.normal(loc=np.asarray(point, dtype=float), scale=effective_rod_length), + dtype=float, + ) + if self.geometry.classify_points(candidate.reshape(1, -1))[0] >= 0: + return candidate + return self.geometry.project_segment_to_domain(np.asarray(point, dtype=float), candidate) + def _apply_point_changes(self, additions: list[FloatArray], removals: list[int]) -> None: """Mutate the point cloud according to the evaluated add/remove plan.""" @@ -201,6 +351,7 @@ def _apply_point_changes(self, additions: list[FloatArray], removals: list[int]) keep[np.asarray(removals, dtype=int)] = False self.points = self.points[keep] self.states = self.states[keep] + self._mark_topology_stale() if additions: additions_arr = _dedupe_points(np.asarray(additions, dtype=float)) @@ -213,6 +364,7 @@ def _apply_point_changes(self, additions: list[FloatArray], removals: list[int]) addition_states, ) ) + self._mark_topology_stale() def _compute_constrained_displacement( self, @@ -253,11 +405,11 @@ def _apply_relaxation_displacements( force_summary: ForceSummary, time_step: float, ) -> float: - """Move dynamic points according to the current force field.""" + """Move dynamic points and return the largest density-scaled relative step.""" new_points = np.array(self.points, copy=True) - max_displacement = 0.0 - max_norm = self.a0 * float(self.params.get("controller_movement_max_freedom", 3.0)) * 0.05 + max_relative_displacement = 0.0 + max_norm = math.inf for index, state in enumerate(self.states): if state not in (STATE_MOBILE, STATE_BOUNDARY): continue @@ -266,15 +418,21 @@ def _apply_relaxation_displacements( displacement = np.array(force_summary.total_forces[index], copy=True) * time_step candidate = self._compute_constrained_displacement(point, int(state), displacement, max_norm) new_points[index] = candidate - max_displacement = max(max_displacement, float(np.linalg.norm(candidate - point))) + density_scale = self.geometry.density_at(candidate) ** (1.0 / max(self.geometry.dim, 1)) + relative_displacement = ( + float(np.linalg.norm(candidate - point)) + * density_scale + / max(self.a0, BOUNDARY_FUZZ) + ) + max_relative_displacement = max(max_relative_displacement, relative_displacement) self.points = new_points - return max_displacement + return max_relative_displacement def _check_convergence(self) -> None: """Update the engine state when a convergence condition is satisfied.""" if ( - self.step >= self.min_equilibrium_steps + self.step > self.min_equilibrium_steps and self.last_max_displacement <= self.tolerated_rel_move ): log.debug( @@ -287,8 +445,7 @@ def _check_convergence(self) -> None: return if ( - self.step >= self.min_equilibrium_steps - and self.last_max_effective_force <= BOUNDARY_FUZZ + self.step >= 2 and self.last_max_effective_force <= BOUNDARY_FUZZ ): log.debug( "Relaxation finished by effective-force equilibrium at step %d (effective_force=%g, rel_move=%g)", @@ -312,17 +469,21 @@ def _step_once(self) -> None: self.a0, self.params, self.step, + simplices=self.current_simplices, ) + self._record_triangulation(force_summary) + self.last_point_density = np.asarray(force_summary.point_density, dtype=float) time_step = self._effective_time_step(force_summary) - max_displacement = self._apply_relaxation_displacements(force_summary, time_step) - self.last_max_displacement = max_displacement / max(self.a0, BOUNDARY_FUZZ) + self.last_max_displacement = self._apply_relaxation_displacements(force_summary, time_step) self.last_max_effective_force = force_summary.max_effective_force self._refresh_boundary_states() + self._update_topology_movement() + self._refresh_topology_if_needed() - if self.step > 0 and self.step % 5 == 0: + if self._should_attempt_point_change(): + self.last_addition_deletion_step = self.step self._attempt_add_delete_points(force_summary) - self._rebuild_mesh() self._check_convergence() def callback_payload(self) -> list[list[Any]]: @@ -343,19 +504,19 @@ def run(self, command: MeshEngineCommand) -> EngineResult: if command != MeshEngineCommand.DO_STEP: raise ValueError(f"Unsupported mesh engine command: {command!r}") - if self.finished or self.step >= self.max_steps: - self._rebuild_mesh() + if self.finished or self._step_limit_reached(): + self._rebuild_mesh(final=True) return MeshEngineStatus.FINISHED_STEP_LIMIT_REACHED, None self.step += 1 self._step_once() if self.finished: - self._rebuild_mesh() + self._rebuild_mesh(final=True) return MeshEngineStatus.FINISHED_FORCE_EQUILIBRIUM_REACHED, None - if self.step >= self.max_steps: - self._rebuild_mesh() + if self._step_limit_reached(): + self._rebuild_mesh(final=True) return MeshEngineStatus.FINISHED_STEP_LIMIT_REACHED, None return MeshEngineStatus.CAN_CONTINUE, self.run @@ -402,5 +563,5 @@ def mesh_bodies_raw( else: engine.run(MeshEngineCommand.DO_STEP) - engine._rebuild_mesh() + engine._rebuild_mesh(final=True) return engine.cached_raw_mesh diff --git a/src/nmesh/mesher/relaxation/finalize.py b/src/nmesh/mesher/relaxation/finalize.py new file mode 100644 index 0000000..3b05a18 --- /dev/null +++ b/src/nmesh/mesher/relaxation/finalize.py @@ -0,0 +1,78 @@ +"""Final mesh cleanup helpers for the relaxation meshing pipeline.""" + +from __future__ import annotations + +import itertools +from typing import Any + +import numpy as np + +from ._constants import BOUNDARY_FUZZ, STATE_BOUNDARY, STATE_FIXED, STATE_MOBILE, STATE_SIMPLE +from ._types import FloatArray +from .geometry import FemGeometry +from .topology import _triangulate_points + + +def snap_final_boundary_points( + points: FloatArray, + states: np.ndarray, + geometry: FemGeometry, + a0: float, + params: dict[str, Any], +) -> tuple[FloatArray, np.ndarray]: + """Snap final points near fixed/boundary neighbors onto implicit boundaries.""" + + if len(points) < geometry.dim + 1: + return points, states + + simplices = _triangulate_points(points, geometry.dim, states) + if len(simplices) == 0: + return points, states + + neighbor_map = _build_neighbor_map(len(points), simplices) + cleaned_points = np.array(points, copy=True) + cleaned_states = np.array(states, copy=True) + acceptable_fuzz = float(params.get("boundary_condition_acceptable_fuzz", BOUNDARY_FUZZ)) + max_steps = int(params.get("boundary_condition_max_nr_correction_steps", 200)) + + for point_index, neighbors in enumerate(neighbor_map): + if not neighbors or not _has_boundary_like_neighbor(cleaned_states, neighbors): + continue + + original = cleaned_points[point_index] + snapped = geometry.project_point_to_boundary_from_inside( + original, + acceptable_fuzz=acceptable_fuzz, + max_steps=max_steps, + ) + if geometry.boundary_distance(snapped) > max(acceptable_fuzz, BOUNDARY_FUZZ): + continue + if geometry.classify_points(snapped.reshape(1, -1))[0] < 0: + continue + local_rod = a0 / (geometry.density_at(original) ** (1.0 / max(geometry.dim, 1))) + if float(np.linalg.norm(snapped - original)) < 0.2 * local_rod: + cleaned_points[point_index] = snapped + if cleaned_states[point_index] in (STATE_MOBILE, STATE_BOUNDARY): + cleaned_states[point_index] = STATE_BOUNDARY + + return cleaned_points, cleaned_states + + +def _build_neighbor_map(point_count: int, simplices: np.ndarray) -> list[list[int]]: + """Build undirected point adjacency from simplices.""" + + neighbors: list[set[int]] = [set() for _ in range(point_count)] + for simplex in simplices: + for left, right in itertools.combinations(simplex.tolist(), 2): + left_index = int(left) + right_index = int(right) + neighbors[left_index].add(right_index) + neighbors[right_index].add(left_index) + return [sorted(group) for group in neighbors] + + +def _has_boundary_like_neighbor(states: np.ndarray, neighbors: list[int]) -> bool: + """Return whether any neighbor should trigger final boundary snapping.""" + + boundary_like_states = {STATE_FIXED, STATE_BOUNDARY, STATE_SIMPLE} + return any(int(states[neighbor]) in boundary_like_states for neighbor in neighbors) diff --git a/src/nmesh/mesher/relaxation/forces.py b/src/nmesh/mesher/relaxation/forces.py index ea728b2..3c6ea2a 100644 --- a/src/nmesh/mesher/relaxation/forces.py +++ b/src/nmesh/mesher/relaxation/forces.py @@ -25,7 +25,13 @@ from ._types import FloatArray from .geometry import FemGeometry from .jit import accumulate_neighbor_forces_default -from .topology import _simplex_measures, _triangulate_points +from .topology import ( + _boundary_state_mask, + _classify_simplices_with_probes, + _simplex_measures, + _simplex_volume_order_ratio, + _triangulate_points, +) @dataclass(frozen=True, slots=True) @@ -50,6 +56,7 @@ class ForceParameters: neigh_force_scale: float irrel_force_scale: float sliver_correction: float + smallest_allowed_volume_ratio: float relaxation_weight: float force_fun: Any boundary_force_fun: Any @@ -210,6 +217,9 @@ def _extract_force_parameters(params: dict[str, Any], step: int) -> ForceParamet neigh_force_scale=float(params.get("controller_neigh_force_scale", 1.0)), irrel_force_scale=float(params.get("controller_irrel_elem_force_scale", 1.0)), sliver_correction=float(params.get("controller_sliver_correction", 1.0)), + smallest_allowed_volume_ratio=float( + params.get("controller_smallest_allowed_volume_ratio", 1.0) + ), relaxation_weight=float(relaxation_weight_fun(step, settling_steps, 0.0, 1.0)), force_fun=params.get("relaxation_force_fun", default_relaxation_force_fun), boundary_force_fun=params.get("boundary_node_force_fun", default_boundary_node_force_fun), @@ -490,14 +500,22 @@ def _compute_simplex_forces( ): return - centroid_regions = geometry.classify_points(np.mean(points[simplices], axis=1)) + relevant_simplices = _classify_relevant_simplices( + points, + states, + geometry, + simplices, + simplex_measures, + dim, + config, + ) suppress_corner_forces = angle_sums < _corner_force_threshold(dim) for simplex_index, simplex in enumerate(simplices): simplex_points = points[simplex] simplex_force = np.zeros_like(simplex_points) volume = float(simplex_measures[simplex_index]) - if centroid_regions[simplex_index] >= 0 and volume > BOUNDARY_FUZZ: - density_here = float(np.mean(point_densities[np.asarray(simplex, dtype=int)])) + if relevant_simplices[simplex_index]: + density_here = geometry.density_at(np.mean(simplex_points, axis=0)) simplex_force = _simplex_force_contribution( simplex_points=simplex_points, dim=dim, @@ -512,9 +530,11 @@ def _compute_simplex_forces( elif config.irrel_force_scale > 0.0: center = np.mean(simplex_points, axis=0) for local_index, point_index in enumerate(simplex.tolist()): - if not _is_dynamic_state(int(states[point_index])): + if int(states[point_index]) != STATE_MOBILE: continue - simplex_force[local_index] = config.irrel_force_scale * (center - simplex_points[local_index]) + simplex_force[local_index] = ( + config.irrel_force_scale * (center - simplex_points[local_index]) + ) for local_index, point_index in enumerate(simplex.tolist()): point_id = int(point_index) @@ -523,6 +543,32 @@ def _compute_simplex_forces( total_forces[point_id] += simplex_force[local_index] +def _classify_relevant_simplices( + points: FloatArray, + states: np.ndarray, + geometry: FemGeometry, + simplices: np.ndarray, + simplex_measures: FloatArray, + dim: int, + config: ForceParameters, +) -> np.ndarray: + """Return the legacy relevant-simplex mask for force calculation.""" + + if len(simplices) == 0: + return np.empty(0, dtype=bool) + + region_ids, probe_consistent = _classify_simplices_with_probes( + points, + simplices, + geometry, + ) + boundary_mask = _boundary_state_mask(points, states, geometry) + all_boundary = np.all(boundary_mask[simplices], axis=1) + boundary_ratio = _simplex_volume_order_ratio(points, simplices, dim) + flat_boundary = all_boundary & (boundary_ratio < config.smallest_allowed_volume_ratio) + return (region_ids >= 0) & probe_consistent & (simplex_measures > BOUNDARY_FUZZ) & ~flat_boundary + + def _finalize_force_summary( total_forces: FloatArray, neighbor_map: list[list[int]], @@ -555,32 +601,35 @@ def _finalize_force_summary( ) boundary_mask = geometry.boundary_mask(points, tolerance=max(0.05 * a0, BOUNDARY_FUZZ)) - full_angle = _full_angle(dim) for point_index, point in enumerate(points): + point_state = int(states[point_index]) + density_here = float(point_densities[point_index]) + effective_rod_length = a0 / (density_here ** (1.0 / max(dim, 1))) + ideal_local_volume = _sphere_volume(0.5 * effective_rod_length, dim) + angle = float(angle_sums[point_index]) point_average_force[point_index] = _average_neighbor_force( neighbor_force_sums[point_index], int(neighbor_force_counts[point_index]), + angle, ) corrected_volume = _corrected_point_volume( simplex_measures=simplex_measures, incident=incident_simplices[point_index], - angle=float(angle_sums[point_index]), - dim=dim, - a0=a0, - is_boundary=bool(boundary_mask[point_index]), - full_angle=full_angle, - ) - density_here = float(point_densities[point_index]) - point_density[point_index] = _point_density_ratio( - density_here=density_here, - corrected_volume=corrected_volume, - a0=a0, + angle=angle, dim=dim, + state=point_state, + has_multiple_immobile_neighbors=_has_multiple_immobile_neighbors( + point_index, + states, + neighbor_map, + ), + ideal_local_volume=ideal_local_volume, ) + point_density[point_index] = ideal_local_volume / max(corrected_volume, DENSITY_EPSILON) point_effective_force[point_index] = _point_effective_force( total_force=total_forces[point_index], - state=int(states[point_index]), + state=point_state, geometry=geometry, point=point, is_boundary=bool(boundary_mask[point_index]), @@ -630,69 +679,74 @@ def _make_force_summary( ) -def _full_angle(dim: int) -> float: - """Return the complete angular measure used for Voronoi-style correction.""" - - if dim == 1: - return 1.5 - if dim == 2: - return 2.0 * math.pi - if dim == 3: - return 4.0 * math.pi - return 1.0 - - -def _average_neighbor_force(force_sum: float, force_count: int) -> float: +def _average_neighbor_force(force_sum: float, force_count: int, angle: float) -> float: """Return the mean absolute neighbor-force magnitude for a point.""" + if angle <= DENSITY_EPSILON: + return 1.0e4 if force_count <= 0: return 0.0 return force_sum / float(force_count) +def _has_multiple_immobile_neighbors( + point_index: int, + states: np.ndarray, + neighbor_map: list[list[int]], +) -> bool: + """Return whether a boundary point is trapped between multiple immobile nodes.""" + + if int(states[point_index]) != STATE_BOUNDARY: + return False + immobile_neighbors = sum( + 1 + for neighbor_index in neighbor_map[point_index] + if int(states[neighbor_index]) in (STATE_FIXED, STATE_SIMPLE) + ) + return immobile_neighbors > 1 + + def _corrected_point_volume( *, simplex_measures: FloatArray, incident: list[int], angle: float, dim: int, - a0: float, - is_boundary: bool, - full_angle: float, + state: int, + has_multiple_immobile_neighbors: bool, + ideal_local_volume: float, ) -> float: """Return the local corrected control volume for one point.""" - if incident: - local_volume = float(np.sum(simplex_measures[np.asarray(incident, dtype=int)])) / float(dim + 1) - else: - local_volume = _sphere_volume(0.5 * a0, dim) + if has_multiple_immobile_neighbors or angle <= DENSITY_EPSILON: + return 1.0e-4 * ideal_local_volume - if angle <= DENSITY_EPSILON: - return 1.0e-4 * _sphere_volume(0.5 * a0, dim) + if not incident: + return 1.0e-4 * ideal_local_volume - correction = ( - full_angle / angle - if dim in (1, 2, 3) - else max(1.0, float(dim + 1) / float(len(incident) or 1)) - ) + incident_measures = simplex_measures[np.asarray(incident, dtype=int)] + if len(incident) == 1: + local_volume = float(incident_measures[0]) / float(max(dim, 1)) + else: + local_volume = float(np.sum(incident_measures)) / float(dim + 1) + + correction = _voronoi_angle_correction(dim, angle, len(incident)) corrected_volume = local_volume * correction - if is_boundary: + if state == STATE_BOUNDARY: corrected_volume *= 1.2 return corrected_volume -def _point_density_ratio( - *, - density_here: float, - corrected_volume: float, - a0: float, - dim: int, -) -> float: - """Return the desired-to-actual local density ratio for one point.""" +def _voronoi_angle_correction(dim: int, angle: float, incident_count: int) -> float: + """Return the legacy angular correction for Voronoi control volumes.""" - effective_rod_length = a0 / (density_here ** (1.0 / max(dim, 1))) - ideal_local_volume = _sphere_volume(0.5 * effective_rod_length, dim) - return ideal_local_volume / max(corrected_volume, DENSITY_EPSILON) + if dim == 1: + return 1.5 + if dim == 2: + return 2.0 * math.pi / angle + if dim == 3: + return 4.0 * math.pi / angle + return math.factorial(dim + 1) / float(max(incident_count, 1)) def _point_effective_force( @@ -725,12 +779,14 @@ def compute_forces( a0: float, params: dict[str, Any], step: int, + simplices: np.ndarray | None = None, ) -> ForceSummary: """Compute neighbor, shape, volume, and irrelevant-element forces.""" point_count = len(points) dim = geometry.dim - simplices = _triangulate_points(points, dim, states) + if simplices is None: + simplices = _triangulate_points(points, dim, states) neighbor_map = _build_neighbor_map(point_count, simplices) config = _extract_force_parameters(params, step) point_densities = np.asarray([geometry.density_at(point) for point in points], dtype=float) diff --git a/src/nmesh/mesher/relaxation/geometry.py b/src/nmesh/mesher/relaxation/geometry.py index da99c21..4e47c46 100644 --- a/src/nmesh/mesher/relaxation/geometry.py +++ b/src/nmesh/mesher/relaxation/geometry.py @@ -144,6 +144,69 @@ def boundary_normal(self, point: FloatArray) -> FloatArray: return best_normal + def boundary_gradient(self, point: FloatArray, body: Body) -> FloatArray: + """Estimate the gradient of one implicit body at a point.""" + + coord = np.asarray(point, dtype=float) + gradient = np.zeros(self.dim, dtype=float) + extent = max(float(np.max(self.bbox_max - self.bbox_min)), 1.0) + epsilon = max(BOUNDARY_FUZZ * 10.0, extent * 1.0e-6) + for axis in range(self.dim): + offset = np.zeros(self.dim, dtype=float) + offset[axis] = epsilon + forward = float(body.evaluate(coord + offset)) + backward = float(body.evaluate(coord - offset)) + gradient[axis] = (forward - backward) / (2.0 * epsilon) + return gradient + + def project_point_to_boundary_from_inside( + self, + point: FloatArray, + *, + acceptable_fuzz: float, + max_steps: int, + ) -> FloatArray: + """Project an interior point onto implicit boundaries using the legacy correction.""" + + coords = np.asarray(point, dtype=float).copy() + bodies = [body for body in self.region_bodies if body is not None] + if not bodies: + return self._project_point_to_box_boundary(coords) + + for _ in range(max(max_steps, 0)): + violated_body: Body | None = None + violated_value = 0.0 + for body in bodies: + value = float(body.evaluate(coords)) + if value > acceptable_fuzz: + violated_body = body + violated_value = value + break + + if violated_body is None: + return coords + + gradient = self.boundary_gradient(coords, violated_body) + gradient_norm_sq = float(np.dot(gradient, gradient)) + scale = 1.0e-6 if gradient_norm_sq <= DENSITY_EPSILON else -violated_value / gradient_norm_sq + coords = coords + scale * gradient + + return coords + + def _project_point_to_box_boundary(self, point: FloatArray) -> FloatArray: + """Project a point to the nearest bounding-box face.""" + + coords = np.asarray(point, dtype=float).copy() + distances_min = np.abs(coords - self.bbox_min) + distances_max = np.abs(coords - self.bbox_max) + min_axis = int(np.argmin(distances_min)) + max_axis = int(np.argmin(distances_max)) + if distances_min[min_axis] <= distances_max[max_axis]: + coords[min_axis] = self.bbox_min[min_axis] + else: + coords[max_axis] = self.bbox_max[max_axis] + return coords + def project_segment_to_domain( self, start: FloatArray, diff --git a/src/nmesh/mesher/relaxation/jit.py b/src/nmesh/mesher/relaxation/jit.py index 57748ea..cb0f202 100644 --- a/src/nmesh/mesher/relaxation/jit.py +++ b/src/nmesh/mesher/relaxation/jit.py @@ -11,7 +11,7 @@ from ._constants import DENSITY_EPSILON, STATE_BOUNDARY, STATE_FIXED, STATE_MOBILE, STATE_SIMPLE -@njit(cache=True) +@njit def _default_relaxation_force(reduced_distance: float) -> float: """Return the legacy default mobile-mobile force law.""" @@ -20,7 +20,7 @@ def _default_relaxation_force(reduced_distance: float) -> float: return 1.0 - reduced_distance -@njit(cache=True) +@njit def _default_boundary_force(reduced_distance: float) -> float: """Return the legacy default boundary interaction force law.""" @@ -31,7 +31,7 @@ def _default_boundary_force(reduced_distance: float) -> float: return 1.0 / reduced_distance - 1.0 -@njit(cache=True) +@njit def accumulate_neighbor_forces_default( points: np.ndarray, states: np.ndarray, diff --git a/src/nmesh/mesher/relaxation/recovery.py b/src/nmesh/mesher/relaxation/recovery.py new file mode 100644 index 0000000..1a7d3e3 --- /dev/null +++ b/src/nmesh/mesher/relaxation/recovery.py @@ -0,0 +1,286 @@ +"""Surface-recovery helpers for the relaxation meshing pipeline.""" + +from __future__ import annotations + +import math +from itertools import combinations + +import numpy as np + +from ._constants import BOUNDARY_FUZZ, DENSITY_EPSILON, STATE_BOUNDARY, STATE_FIXED, STATE_MOBILE, STATE_SIMPLE +from ._types import FloatArray +from .geometry import FemGeometry +from .topology import _simplex_measures + + +def mirror_surface_recovery_points( + points: FloatArray, + states: np.ndarray, + geometry: FemGeometry, + a0: float, + simplices: np.ndarray, +) -> FloatArray: + """Return legacy mirror/prevention points for poor boundary simplices.""" + + dim = geometry.dim + if dim not in (2, 3) or len(simplices) == 0: + return np.empty((0, dim), dtype=float) + + simplex_measures = _simplex_measures(points, simplices, dim) + seen_surfaces: set[tuple[tuple[float, ...], ...]] = set() + additions: list[FloatArray] = [] + + for simplex, simplex_volume in zip(simplices, simplex_measures): + simplex_points = np.asarray(points[simplex], dtype=float) + boundary_points, interior_points = _split_boundary_and_interior_points( + simplex, + simplex_points, + states, + ) + if dim == 2: + additions.extend( + _mirror_points_2d( + boundary_points, + interior_points, + simplex_points, + float(simplex_volume), + geometry, + a0, + seen_surfaces, + ) + ) + else: + additions.extend( + _prevention_points_3d( + boundary_points, + simplex_points, + float(simplex_volume), + geometry, + a0, + seen_surfaces, + ) + ) + + if not additions: + return np.empty((0, dim), dtype=float) + return _dedupe_recovery_points(np.asarray(additions, dtype=float)) + + +def _split_boundary_and_interior_points( + simplex: np.ndarray, + simplex_points: FloatArray, + states: np.ndarray, +) -> tuple[list[FloatArray], list[FloatArray]]: + """Split simplex vertices into boundary-like and mobile interior vertices.""" + + boundary_points: list[FloatArray] = [] + interior_points: list[FloatArray] = [] + for local_index, point_index in enumerate(simplex.tolist()): + state = int(states[int(point_index)]) + point = simplex_points[local_index] + if state in (STATE_BOUNDARY, STATE_FIXED, STATE_SIMPLE): + boundary_points.append(point) + elif state == STATE_MOBILE: + interior_points.append(point) + return boundary_points, interior_points + + +def _mirror_points_2d( + boundary_points: list[FloatArray], + interior_points: list[FloatArray], + simplex_points: FloatArray, + simplex_volume: float, + geometry: FemGeometry, + a0: float, + seen_surfaces: set[tuple[tuple[float, ...], ...]], +) -> list[FloatArray]: + """Return mirrored points for the legacy 2D boundary-simplex cases.""" + + additions: list[FloatArray] = [] + if len(boundary_points) == 3: + shifted = list(boundary_points) + for _ in range(3): + vertex_point = shifted[0] + surface_points = np.asarray(shifted[1:], dtype=float) + if _register_surface(surface_points, seen_surfaces): + shifted = shifted[1:] + shifted[:1] + continue + additions.extend( + _node_action_2d(surface_points, vertex_point, simplex_points, simplex_volume, geometry, a0) + ) + shifted = shifted[1:] + shifted[:1] + return additions + + if len(boundary_points) == 2 and len(interior_points) == 1: + surface_points = np.asarray(boundary_points, dtype=float) + if not _register_surface(surface_points, seen_surfaces): + additions.extend( + _node_action_2d( + surface_points, + interior_points[0], + simplex_points, + simplex_volume, + geometry, + a0, + ) + ) + return additions + + +def _node_action_2d( + surface_points: FloatArray, + vertex_point: FloatArray, + simplex_points: FloatArray, + simplex_volume: float, + geometry: FemGeometry, + a0: float, +) -> list[FloatArray]: + """Mirror one 2D point when the surface angle and volume tests request it.""" + + angle = _solid_angle_2d(surface_points, vertex_point) + if angle > 0.55 * math.pi and _volume_ratio(simplex_points, simplex_volume, geometry, a0) > 3.0e-2: + return [_mirror_point(surface_points, vertex_point, 2)] + return [] + + +def _prevention_points_3d( + boundary_points: list[FloatArray], + simplex_points: FloatArray, + simplex_volume: float, + geometry: FemGeometry, + a0: float, + seen_surfaces: set[tuple[tuple[float, ...], ...]], +) -> list[FloatArray]: + """Return midpoint prevention points for legacy 3D boundary edges.""" + + if len(boundary_points) < 2: + return [] + + additions: list[FloatArray] = [] + for left, right in combinations(boundary_points, 2): + surface_points = np.asarray([left, right], dtype=float) + already_seen = _register_surface(surface_points, seen_surfaces) + if already_seen: + continue + middle_point = 0.5 * (surface_points[0] + surface_points[1]) + middle_density = geometry.density_at(middle_point) + middle_rod_length = a0 / (middle_density ** (1.0 / max(geometry.dim, 1))) + distance = float(np.linalg.norm(surface_points[0] - surface_points[1])) + if ( + distance > 1.2 * middle_rod_length + and _volume_ratio(simplex_points, simplex_volume, geometry, a0) > 3.0e-2 + and _is_outer_region_constrained(geometry, middle_point) + ): + additions.append(middle_point) + return additions + + +def _surface_key(surface_points: FloatArray) -> tuple[tuple[float, ...], ...]: + """Return a stable set-like key for a boundary surface.""" + + return tuple(sorted(tuple(np.round(point, decimals=10).tolist()) for point in surface_points)) + + +def _register_surface( + surface_points: FloatArray, + seen_surfaces: set[tuple[tuple[float, ...], ...]], +) -> bool: + """Register a surface and return whether it had already been processed.""" + + key = _surface_key(surface_points) + if key in seen_surfaces: + return True + seen_surfaces.add(key) + return False + + +def _solid_angle_2d(surface_points: FloatArray, vertex_point: FloatArray) -> float: + """Return the angle subtended by two surface points from a vertex.""" + + first = surface_points[0] - vertex_point + second = surface_points[1] - vertex_point + denom = float(np.linalg.norm(first) * np.linalg.norm(second)) + if denom <= DENSITY_EPSILON: + return 0.0 + cosine = float(np.clip(np.dot(first, second) / denom, -1.0, 1.0)) + return float(math.acos(cosine)) + + +def _mirror_point(surface_points: FloatArray, point: FloatArray, dim: int) -> FloatArray: + """Mirror a point across a line/plane defined by boundary points.""" + + if len(surface_points) == dim: + origin = surface_points[0] + directions = surface_points[1:] - origin + _, _, vh = np.linalg.svd(directions, full_matrices=True) + normal = np.asarray(vh[-1], dtype=float) + normal_norm_sq = float(np.dot(normal, normal)) + if normal_norm_sq <= DENSITY_EPSILON: + return np.asarray(point, dtype=float) + factor = float(np.dot(point - origin, normal)) / normal_norm_sq + projection = point - factor * normal + return np.asarray(2.0 * projection - point, dtype=float) + + if len(surface_points) == dim - 1 and dim == 3: + start = surface_points[0] + segment = surface_points[1] - start + denom = float(np.dot(segment, segment)) + if denom <= DENSITY_EPSILON: + return np.asarray(point, dtype=float) + fraction = float(np.dot(point - start, segment)) / denom + projection = start + fraction * segment + return np.asarray(2.0 * projection - point, dtype=float) + + return np.asarray(point, dtype=float) + + +def _volume_ratio( + simplex_points: FloatArray, + simplex_volume: float, + geometry: FemGeometry, + a0: float, +) -> float: + """Return the legacy real-to-ideal simplex volume ratio.""" + + if simplex_volume <= DENSITY_EPSILON: + return 0.0 + center = np.mean(simplex_points, axis=0) + density_here = geometry.density_at(center) + rod_scaled = a0 / (density_here ** (1.0 / max(geometry.dim, 1))) + ideal_volume = _regular_simplex_volume(rod_scaled, geometry.dim) + return abs(simplex_volume / max(ideal_volume, DENSITY_EPSILON)) + + +def _regular_simplex_volume(edge_length: float, dim: int) -> float: + """Return the volume of a regular simplex with the supplied edge length.""" + + numerator = edge_length**dim * math.sqrt(dim + 1.0) + denominator = math.factorial(dim) * math.sqrt(2.0**dim) + return numerator / denominator + + +def _is_outer_region_constrained(geometry: FemGeometry, point: FloatArray) -> bool: + """Return whether a point satisfies the legacy outer-region boundary check.""" + + if not geometry.points_in_bbox(np.asarray(point, dtype=float).reshape(1, -1))[0]: + return False + for body in geometry.region_bodies: + if body is None: + continue + if float(body.evaluate(point)) >= BOUNDARY_FUZZ: + return False + return True + + +def _dedupe_recovery_points(points: FloatArray) -> FloatArray: + """Deduplicate generated recovery points while preserving order.""" + + keep_indices: list[int] = [] + seen: set[tuple[float, ...]] = set() + for index, point in enumerate(points): + key = tuple(np.round(point, decimals=10).tolist()) + if key in seen: + continue + seen.add(key) + keep_indices.append(index) + return points[np.asarray(keep_indices, dtype=int)] diff --git a/src/nmesh/mesher/relaxation/seeding.py b/src/nmesh/mesher/relaxation/seeding.py index 0c5a5b1..b747c0e 100644 --- a/src/nmesh/mesher/relaxation/seeding.py +++ b/src/nmesh/mesher/relaxation/seeding.py @@ -2,13 +2,13 @@ from __future__ import annotations +import math from typing import Any import numpy as np from ._constants import ( BOUNDARY_FUZZ, - DEFAULT_RNG_SEED, DENSITY_EPSILON, STATE_BOUNDARY, STATE_FIXED, @@ -95,45 +95,129 @@ def filter_points(points: FloatArray) -> FloatArray: ) -def _grid_candidate_points(geometry: FemGeometry, a0: float) -> FloatArray: - """Generate deterministic candidate seed points over the bounding box.""" +def _box_volume(geometry: FemGeometry) -> float: + """Return the volume of the meshing bounding box.""" - extents = np.maximum(geometry.bbox_max - geometry.bbox_min, a0) - sample_points = [geometry.bbox_min, geometry.bbox_max, 0.5 * (geometry.bbox_min + geometry.bbox_max)] - max_density = max(geometry.density_at(point) for point in sample_points) - density_scale = min(max_density ** (1.0 / max(geometry.dim, 1)), 2.0) - step = max(a0 / density_scale, a0 * 0.5) + return float(np.prod(np.maximum(geometry.bbox_max - geometry.bbox_min, 0.0))) - counts = np.maximum(np.floor(extents / step).astype(int) + 1, 2) - candidate_count = int(np.prod(counts)) - if candidate_count > 15_000: - approximate = min(8_000, max(500, int(np.prod(extents / max(a0, DENSITY_EPSILON)) * 4))) - rng = np.random.default_rng(DEFAULT_RNG_SEED) - return rng.uniform(geometry.bbox_min, geometry.bbox_max, size=(approximate, geometry.dim)) - axes = [ - np.linspace(geometry.bbox_min[axis], geometry.bbox_max[axis], int(counts[axis])) - for axis in range(geometry.dim) - ] - grid = np.meshgrid(*axes, indexing="ij") - return np.stack(grid, axis=-1).reshape(-1, geometry.dim) +def _random_point_in_box(geometry: FemGeometry, rng: np.random.Generator) -> FloatArray: + """Sample one uniformly random point from the bounding box.""" + return geometry.bbox_min + rng.random(geometry.dim) * (geometry.bbox_max - geometry.bbox_min) -def _keep_point_for_density(point: FloatArray, density_here: float, rng: np.random.Generator) -> bool: - """Return whether a candidate point survives density-based thinning.""" - _ = point - if density_here >= 1.0: - return True - return bool(rng.random() <= density_here) +def _sampling_density(geometry: FemGeometry, point: FloatArray) -> float: + """Evaluate density only for points belonging to the meshed domain.""" + + if geometry.classify_points(np.asarray(point, dtype=float).reshape(1, -1))[0] < 0: + return 0.0 + return geometry.density_at(point) + + +def _estimate_density_max_and_average( + geometry: FemGeometry, + nr_probes: int, + rng: np.random.Generator, + *, + conservative_factor: float = 1.12, +) -> tuple[float, float]: + """Estimate the legacy random-sampling max and average density values.""" + + if nr_probes <= 0: + return 0.0, 0.0 + + first_value = _sampling_density(geometry, _random_point_in_box(geometry, rng)) + max_seen = first_value + # The legacy estimator intentionally seeded the max with the first probe + # while leaving it out of the average accumulator. + density_sum = 0.0 + for _ in range(1, nr_probes): + value = _sampling_density(geometry, _random_point_in_box(geometry, rng)) + max_seen = max(max_seen, value) + density_sum += value + return max_seen * conservative_factor, density_sum / float(nr_probes) + + +def _sphere_volume(dim: int) -> float: + """Return the d-dimensional unit sphere volume.""" + + return (math.pi ** (0.5 * dim)) / math.gamma(1.0 + 0.5 * dim) + + +def _sphere_packing_ratio_lattice_type_d(dim: int) -> float: + """Return the D-lattice sphere-packing ratio used by the legacy seeder.""" + + if dim <= 1: + return 1.0 + + lattice_vectors = np.zeros((dim, dim), dtype=float) + for row in range(dim): + if row == dim - 1: + lattice_vectors[row, max(dim - 2, 0) :] = 1.0 + else: + lattice_vectors[row, row] = 1.0 + lattice_vectors[row, row + 1] = -1.0 + + lattice_cell_volume = abs(float(np.linalg.det(lattice_vectors))) + if lattice_cell_volume <= DENSITY_EPSILON: + return 1.0 + sphere_radius = 0.5 * math.sqrt(2.0) + sphere_volume = (sphere_radius**dim) * _sphere_volume(dim) + return sphere_volume / lattice_cell_volume + + +def _estimate_initial_point_count( + geometry: FemGeometry, + a0: float, + fixed_points: FloatArray, + average_density: float, +) -> int: + """Estimate the number of random initial nodes using the OCaml formula.""" + + _ = fixed_points + node_volume = ( + _sphere_volume(geometry.dim) + * ((a0 * 0.5 * 0.7) ** geometry.dim) + / max(_sphere_packing_ratio_lattice_type_d(geometry.dim), DENSITY_EPSILON) + ) + estimated_integrated_density = average_density * _box_volume(geometry) + estimated_nodes = estimated_integrated_density / max(node_volume, DENSITY_EPSILON) + return min(10_000, max(geometry.dim + 1 + 5, int(estimated_nodes))) + + +def _distribute_points_randomly( + geometry: FemGeometry, + nr_points: int, + max_density: float, + rng: np.random.Generator, +) -> FloatArray: + """Draw initial points with rejection sampling against the density field.""" + + if nr_points <= 0 or max_density <= DENSITY_EPSILON: + return np.empty((0, geometry.dim), dtype=float) + + # OCaml's Array.make evaluates the initial random point once before the + # rejection loop; keep the same RNG consumption for deterministic parity. + _ = _random_point_in_box(geometry, rng) + result = np.empty((nr_points, geometry.dim), dtype=float) + accepted = 0 + while accepted < nr_points: + point = _random_point_in_box(geometry, rng) + if rng.random() * max_density > _sampling_density(geometry, point): + continue + result[accepted] = point + accepted += 1 + return result def _classify_dynamic_states(geometry: FemGeometry, points: FloatArray, a0: float) -> np.ndarray: """Return mobile or boundary states for the supplied movable points.""" + _ = a0 if len(points) == 0: return np.empty(0, dtype=int) - mask = geometry.boundary_mask(points, tolerance=max(0.05 * a0, BOUNDARY_FUZZ)) + mask = geometry.boundary_mask(points, tolerance=BOUNDARY_FUZZ) states = np.full(len(points), STATE_MOBILE, dtype=int) states[mask] = STATE_BOUNDARY return states @@ -155,12 +239,14 @@ def _select_generated_points( mobile_points: FloatArray, simply_points: FloatArray, rng: np.random.Generator, + params: dict[str, Any], ) -> FloatArray: - """Generate and thin candidate points that are not already seeded.""" + """Generate density-weighted random seed points like the legacy mesher.""" - candidates = _grid_candidate_points(geometry, a0) - region_ids = geometry.classify_points(candidates) - candidates = candidates[region_ids >= 0] + nr_probes = int(params.get("nr_probes_for_determining_volume", 100_000)) + max_density, average_density = _estimate_density_max_and_average(geometry, nr_probes, rng) + nr_points = _estimate_initial_point_count(geometry, a0, fixed_points, average_density) + candidates = _distribute_points_randomly(geometry, nr_points, max_density, rng) if len(candidates) == 0: return np.empty((0, geometry.dim), dtype=float) @@ -173,8 +259,7 @@ def _select_generated_points( key = _point_key(point) if key in candidate_keys: continue - density_here = geometry.density_at(point) - if not _keep_point_for_density(point, density_here, rng): + if geometry.classify_points(np.asarray(point, dtype=float).reshape(1, -1))[0] < 0: continue selected_candidates.append(point) candidate_keys.add(key) @@ -194,6 +279,48 @@ def _collect_hint_points(geometry: FemGeometry) -> FloatArray: return _filter_relevant_points(geometry, hint_block) +def _periodic_outer_box_points( + geometry: FemGeometry, + a0: float, + periodic: list[float] | list[bool], +) -> FloatArray: + """Create paired fixed seed points on periodic outer-box faces.""" + + periodic_flags = [bool(value) for value in periodic] + if not any(periodic_flags): + return np.empty((0, geometry.dim), dtype=float) + + points: list[FloatArray] = [] + spacing = max(a0, DENSITY_EPSILON) + for periodic_axis, enabled in enumerate(periodic_flags): + if not enabled: + continue + + other_axes = [axis for axis in range(geometry.dim) if axis != periodic_axis] + face_axes: list[FloatArray] = [] + for axis in other_axes: + extent = float(geometry.bbox_max[axis] - geometry.bbox_min[axis]) + count = max(2, int(math.floor(extent / spacing)) + 1) + face_axes.append(np.linspace(geometry.bbox_min[axis], geometry.bbox_max[axis], count)) + + coordinate_rows = ( + np.zeros((1, 0), dtype=float) + if not other_axes + else np.array(np.meshgrid(*face_axes, indexing="ij")).T.reshape(-1, len(other_axes)) + ) + for coordinates in coordinate_rows: + for side in (geometry.bbox_min[periodic_axis], geometry.bbox_max[periodic_axis]): + point = np.zeros(geometry.dim, dtype=float) + point[periodic_axis] = side + for axis, value in zip(other_axes, coordinates): + point[axis] = value + points.append(point) + + if not points: + return np.empty((0, geometry.dim), dtype=float) + return _filter_relevant_points(geometry, _dedupe_points(np.asarray(points, dtype=float))) + + def _apply_periodic_fixed_states( states: np.ndarray, all_points: FloatArray, @@ -208,7 +335,8 @@ def _apply_periodic_fixed_states( periodic_flags = [bool(value) for value in periodic] periodic_mask = np.zeros(len(all_points), dtype=bool) - tolerance = max(0.05 * a0, BOUNDARY_FUZZ) + _ = a0 + tolerance = BOUNDARY_FUZZ for axis, enabled in enumerate(periodic_flags): if not enabled: continue @@ -229,15 +357,23 @@ def _prepare_initial_points( simply_points: FloatArray, periodic: list[float] | list[bool], rng: np.random.Generator, + params: dict[str, Any] | None = None, ) -> tuple[FloatArray, np.ndarray]: """Prepare the initial point cloud and point-state array for relaxation.""" + params = {} if params is None else params fixed_points, mobile_points, simply_points = _dedupe_fixed_mobile( fixed_points, mobile_points, simply_points ) fixed_points = _filter_relevant_points(geometry, fixed_points) mobile_points = _filter_relevant_points(geometry, mobile_points) simply_points = _filter_relevant_points(geometry, simply_points) + + if len(simply_points) > 0: + states = np.full(len(simply_points), STATE_SIMPLE, dtype=int) + _apply_periodic_fixed_states(states, simply_points, geometry, periodic, a0) + return simply_points, states + generated_points = _select_generated_points( geometry, a0, @@ -245,8 +381,10 @@ def _prepare_initial_points( mobile_points, simply_points, rng, - ) + params, + ) if len(mobile_points) == 0 else np.empty((0, geometry.dim), dtype=float) hint_block = _collect_hint_points(geometry) + periodic_block = _periodic_outer_box_points(geometry, a0, periodic) all_points = np.vstack( ( @@ -254,6 +392,7 @@ def _prepare_initial_points( simply_points, mobile_points, hint_block, + periodic_block, generated_points, ) ) @@ -264,6 +403,7 @@ def _prepare_initial_points( np.full(len(simply_points), STATE_SIMPLE, dtype=int), _classify_dynamic_states(geometry, mobile_points, a0), np.full(len(hint_block), STATE_FIXED, dtype=int), + np.full(len(periodic_block), STATE_FIXED, dtype=int), _classify_dynamic_states(geometry, generated_points, a0), ) ) diff --git a/src/nmesh/mesher/relaxation/topology.py b/src/nmesh/mesher/relaxation/topology.py index 2b42674..d02c07d 100644 --- a/src/nmesh/mesher/relaxation/topology.py +++ b/src/nmesh/mesher/relaxation/topology.py @@ -173,6 +173,24 @@ def _triangulate_points(points: FloatArray, dim: int, states: np.ndarray | None return np.empty((0, dim + 1), dtype=int) +def _orient_simplices_positive(points: FloatArray, simplices: np.ndarray, dim: int) -> np.ndarray: + """Return simplices ordered with positive legacy simplex orientation.""" + + if dim <= 1 or len(simplices) == 0: + return simplices + + oriented = np.array(simplices, copy=True) + matrices = points[oriented[:, 1:]] - points[oriented[:, [0]]] + determinants = np.linalg.det(matrices) + negative = np.flatnonzero(determinants < 0.0) + if len(negative) > 0: + oriented[negative, 0], oriented[negative, 1] = ( + oriented[negative, 1].copy(), + oriented[negative, 0].copy(), + ) + return oriented + + def assemble_raw_mesh( points: FloatArray, geometry: FemGeometry, @@ -186,21 +204,22 @@ def assemble_raw_mesh( coords = np.asarray(points, dtype=float) dim = geometry.dim simplices = _triangulate_points(coords, dim, states) + simplices = _orient_simplices_positive(coords, simplices, dim) if len(simplices) > 0: region_ids, probe_consistent = _classify_simplices_with_probes(coords, simplices, geometry) measures = _simplex_measures(coords, simplices, dim) - boundary_mask = geometry.boundary_mask( + boundary_mask = _boundary_state_mask( coords, - tolerance=max(BOUNDARY_FUZZ * 10.0, 1.0e-5), + states, + geometry, ) all_boundary = np.all(boundary_mask[simplices], axis=1) boundary_ratio = _simplex_volume_order_ratio(coords, simplices, dim) - normalized_ratio = boundary_ratio / max(_regular_boundary_ratio(dim), BOUNDARY_FUZZ) smallest_allowed_ratio = float( (params or {}).get("controller_smallest_allowed_volume_ratio", 1.0) ) - flat_boundary = all_boundary & (normalized_ratio < smallest_allowed_ratio) + flat_boundary = all_boundary & (boundary_ratio < smallest_allowed_ratio) keep = probe_consistent & (measures > BOUNDARY_FUZZ) & ~flat_boundary simplices = simplices[keep] region_ids = region_ids[keep] @@ -235,6 +254,21 @@ def assemble_raw_mesh( ) +def _boundary_state_mask( + coords: FloatArray, + states: np.ndarray | None, + geometry: FemGeometry, +) -> np.ndarray: + """Return the point mask used for legacy flat-boundary simplex rejection.""" + + if states is not None and len(states) == len(coords): + return np.asarray(states, dtype=int) == STATE_BOUNDARY + return geometry.boundary_mask( + coords, + tolerance=max(BOUNDARY_FUZZ * 10.0, 1.0e-5), + ) + + def _callback_mesh_info(raw_mesh: RawMesh) -> list[list[Any]]: """Build the legacy-style callback payload expected by old consumers.""" diff --git a/tests/nmesh/test_integration.py b/tests/nmesh/test_integration.py index 48eb539..eed9dea 100644 --- a/tests/nmesh/test_integration.py +++ b/tests/nmesh/test_integration.py @@ -8,6 +8,8 @@ def test_deterministic_meshing_uses_fixed_rng_seed(): bounding_box=[[0.0, 0.0], [1.0, 1.0]], objects=[nmesh.Box([0.1, 0.1], [0.9, 0.9])], a0=0.35, + max_steps=20, + nr_probes_for_determining_volume=500, ) mesh_a = nmesh.Mesh(**kwargs) @@ -25,6 +27,8 @@ def test_periodic_scenario_completes_with_valid_topology(): mesh_bounding_box=True, periodic=[True, False], a0=0.4, + max_steps=20, + nr_probes_for_determining_volume=500, ) assert len(mesh.points) > 0 @@ -34,6 +38,23 @@ def test_periodic_scenario_completes_with_valid_topology(): assert all(region >= 0 for region in mesh.regions) +def test_multi_axis_periodic_mesh_groups_corners_and_edges(): + mesh = nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[nmesh.Box([0.0, 0.0], [1.0, 1.0])], + mesh_bounding_box=True, + periodic=[True, True], + a0=0.4, + max_steps=20, + nr_probes_for_determining_volume=600, + ) + + group_sizes = sorted(len(group) for group in mesh.periodic_point_indices) + assert len(mesh.simplices) > 0 + assert 4 in group_sizes + assert group_sizes.count(2) >= 2 + + def test_hint_driven_scenario_completes_with_valid_topology(): seed_mesh = nmesh.mesh_from_points_and_simplices( points=[[0.2, 0.2], [0.8, 0.2], [0.5, 0.75]], @@ -47,6 +68,8 @@ def test_hint_driven_scenario_completes_with_valid_topology(): mesh_bounding_box=True, hints=[(seed_mesh, nmesh.Box([0.1, 0.1], [0.9, 0.9]))], a0=0.4, + max_steps=20, + nr_probes_for_determining_volume=500, ) points = np.asarray(mesh.points, dtype=float) @@ -56,3 +79,44 @@ def test_hint_driven_scenario_completes_with_valid_topology(): assert np.any(np.all(np.isclose(points, [0.2, 0.2]), axis=1)) assert np.any(np.all(np.isclose(points, [0.8, 0.2]), axis=1)) assert np.any(np.all(np.isclose(points, [0.5, 0.75]), axis=1)) + + +def test_adjacent_regions_keep_valid_piece_topology(): + mesh = nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[ + nmesh.Box([0.0, 0.0], [0.5, 1.0]), + nmesh.Box([0.5, 0.0], [1.0, 1.0]), + ], + a0=0.5, + max_steps=20, + nr_probes_for_determining_volume=500, + ) + + assert len(mesh.points) > 0 + assert len(mesh.simplices) > 0 + assert sorted(set(mesh.regions)) == [1, 2] + assert len(mesh.surfaces) > 0 + + +def test_concave_difference_domain_keeps_hole_empty(): + hole = nmesh.Box([0.35, 0.35], [0.65, 0.65]) + shell = nmesh.difference( + nmesh.Box([0.0, 0.0], [1.0, 1.0]), + [hole], + ) + + mesh = nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[shell], + a0=0.25, + max_steps=20, + nr_probes_for_determining_volume=1000, + ) + + points = np.asarray(mesh.points, dtype=float) + centroids = np.mean(points[np.asarray(mesh.simplices, dtype=int)], axis=1) + inside_hole = np.asarray(hole.obj.evaluate(centroids), dtype=float) > 0.0 + assert len(mesh.simplices) > 0 + assert not np.any(inside_hole) + assert len(mesh.surfaces) >= 8 diff --git a/tests/nmesh/test_meshing_engine.py b/tests/nmesh/test_meshing_engine.py index 1df9785..8986a8a 100644 --- a/tests/nmesh/test_meshing_engine.py +++ b/tests/nmesh/test_meshing_engine.py @@ -7,20 +7,33 @@ _compile_density_function, assemble_raw_mesh, fem_geometry_from_bodies, + RelaxationEngine, ) -from nmesh.mesher.relaxation._constants import STATE_BOUNDARY, STATE_MOBILE -from nmesh.mesher.relaxation.forces import compute_forces -from nmesh.mesher.relaxation.seeding import _classify_dynamic_states +from nmesh.mesher.relaxation._constants import BOUNDARY_FUZZ, STATE_BOUNDARY, STATE_MOBILE, STATE_SIMPLE +from nmesh.mesher.relaxation._constants import STATE_FIXED +from nmesh.mesher.relaxation.finalize import snap_final_boundary_points +from nmesh.mesher.relaxation.forces import ( + _classify_relevant_simplices, + _extract_force_parameters, + compute_forces, +) +from nmesh.mesher.relaxation.recovery import mirror_surface_recovery_points +from nmesh.mesher.relaxation.seeding import _classify_dynamic_states, _prepare_initial_points +from nmesh.mesher.relaxation.topology import _orient_simplices_positive +from nmesh.mesher.periodic import build_periodic_groups class TestMeshingEngine(unittest.TestCase): - def test_boundary_state_detection_marks_surface_points(self): - geometry = fem_geometry_from_bodies( + def _box_geometry(self): + return fem_geometry_from_bodies( (np.asarray([0.0, 0.0]), np.asarray([1.0, 1.0])), [Box([0.0, 0.0], [1.0, 1.0]).obj], [], ) + def test_boundary_state_detection_marks_surface_points(self): + geometry = self._box_geometry() + points = np.asarray( [ [0.0, 0.5], @@ -34,11 +47,7 @@ def test_boundary_state_detection_marks_surface_points(self): self.assertEqual(states.tolist(), [STATE_BOUNDARY, STATE_MOBILE]) def test_force_summary_suppresses_shape_force_at_boundary_corners(self): - geometry = fem_geometry_from_bodies( - (np.asarray([0.0, 0.0]), np.asarray([1.0, 1.0])), - [Box([0.0, 0.0], [1.0, 1.0]).obj], - [], - ) + geometry = self._box_geometry() points = np.asarray( [ @@ -73,11 +82,7 @@ def test_force_summary_suppresses_shape_force_at_boundary_corners(self): self.assertGreater(summary.max_effective_force, 0.0) def test_boundary_effective_force_uses_tangential_projection(self): - geometry = fem_geometry_from_bodies( - (np.asarray([0.0, 0.0]), np.asarray([1.0, 1.0])), - [Box([0.0, 0.0], [1.0, 1.0]).obj], - [], - ) + geometry = self._box_geometry() points = np.asarray( [ @@ -108,12 +113,57 @@ def test_boundary_effective_force_uses_tangential_projection(self): self.assertAlmostEqual(summary.total_forces[0][1], 0.0, places=7) self.assertAlmostEqual(summary.point_effective_force[0], 0.0, places=7) - def test_geometry_projects_outside_point_back_to_boundary(self): - geometry = fem_geometry_from_bodies( - (np.asarray([0.0, 0.0]), np.asarray([1.0, 1.0])), - [Box([0.0, 0.0], [1.0, 1.0]).obj], - [], + def test_force_relevant_simplex_classifier_rejects_boundary_crossing_probe(self): + geometry = self._box_geometry() + points = np.asarray( + [ + [0.5, 0.5], + [1.2, 0.5], + [0.5, 0.8], + ], + dtype=float, ) + simplices = np.asarray([[0, 1, 2]], dtype=int) + states = np.full(len(points), STATE_MOBILE, dtype=int) + measures = np.asarray([0.105], dtype=float) + + mask = _classify_relevant_simplices( + points, + states, + geometry, + simplices, + measures, + 2, + _extract_force_parameters({}, step=10), + ) + + self.assertFalse(bool(mask[0])) + + def test_surface_recovery_mirrors_poor_2d_boundary_simplex(self): + geometry = self._box_geometry() + points = np.asarray( + [ + [0.0, 0.0], + [1.0, 0.0], + [0.5, 0.2], + ], + dtype=float, + ) + states = np.asarray([STATE_BOUNDARY, STATE_BOUNDARY, STATE_MOBILE], dtype=int) + + recovery_points = mirror_surface_recovery_points( + points, + states, + geometry, + a0=0.5, + simplices=np.asarray([[0, 1, 2]], dtype=int), + ) + + self.assertEqual(len(recovery_points), 1) + self.assertTrue(np.allclose(recovery_points[0], [0.5, -0.2])) + + def test_geometry_projects_outside_point_back_to_boundary(self): + geometry = self._box_geometry() projected = geometry.project_segment_to_domain( np.asarray([0.2, 0.5], dtype=float), @@ -123,6 +173,63 @@ def test_geometry_projects_outside_point_back_to_boundary(self): self.assertGreaterEqual(geometry.classify_points(projected.reshape(1, -1))[0], 0) self.assertAlmostEqual(projected[0], 0.0, places=4) + def test_final_cleanup_snaps_near_fixed_neighbors_to_boundary(self): + geometry = self._box_geometry() + points = np.asarray( + [ + [0.0, 0.0], + [0.0, 1.0], + [0.05, 0.5], + ], + dtype=float, + ) + states = np.asarray([STATE_FIXED, STATE_FIXED, STATE_MOBILE], dtype=int) + + snapped_points, snapped_states = snap_final_boundary_points( + points, + states, + geometry, + a0=0.5, + params={ + "boundary_condition_acceptable_fuzz": 1.0e-6, + "boundary_condition_max_nr_correction_steps": 200, + }, + ) + + self.assertAlmostEqual(snapped_points[2][0], 0.0, places=5) + self.assertEqual(int(snapped_states[2]), STATE_BOUNDARY) + + def test_final_cleanup_removes_outside_dynamic_points(self): + geometry = self._box_geometry() + engine = RelaxationEngine( + geometry, + {"parameters": {}}, + 0.5, + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + [False, False], + ) + engine.points = np.asarray( + [ + [0.2, 0.2], + [0.8, 0.2], + [0.2, 0.8], + [1.2, 0.5], + ], + dtype=float, + ) + engine.states = np.asarray( + [STATE_MOBILE, STATE_MOBILE, STATE_MOBILE, STATE_MOBILE], + dtype=int, + ) + engine.last_point_density = np.ones(len(engine.points), dtype=float) + + points, states = engine._final_output_points_and_states() + + self.assertEqual(len(points), 3) + self.assertEqual(states.tolist(), [STATE_MOBILE, STATE_MOBILE, STATE_MOBILE]) + def test_assemble_raw_mesh_filters_flat_boundary_sliver(self): geometry = fem_geometry_from_bodies( (np.asarray([0.0, 0.0]), np.asarray([1.0, 0.01])), @@ -170,6 +277,65 @@ def test_assemble_raw_mesh_keeps_regular_boundary_simplex(self): self.assertEqual(len(raw_mesh.simplices), 1) + def test_simplex_orientation_is_normalized_positive(self): + points = np.asarray( + [ + [0.0, 0.0], + [1.0, 0.0], + [0.0, 1.0], + ], + dtype=float, + ) + simplices = np.asarray([[0, 2, 1]], dtype=int) + + oriented = _orient_simplices_positive(points, simplices, dim=2) + determinant = np.linalg.det(points[oriented[0, 1:]] - points[oriented[0, [0]]]) + + self.assertGreater(determinant, 0.0) + + def test_periodic_groups_merge_multi_axis_corners(self): + points = np.asarray( + [ + [0.0, 0.0], + [0.0, 1.0], + [1.0, 0.0], + [1.0, 1.0], + [0.5, 0.0], + [0.5, 1.0], + ], + dtype=float, + ) + + groups = build_periodic_groups( + points, + np.asarray([0.0, 0.0]), + np.asarray([1.0, 1.0]), + [True, True], + tolerance=1.0e-6, + ) + + self.assertIn([0, 1, 2, 3], groups) + self.assertIn([4, 5], groups) + + def test_periodic_groups_do_not_round_non_periodic_coordinates(self): + points = np.asarray( + [ + [0.0, 0.5], + [1.0, 0.5 + 1.0e-9], + ], + dtype=float, + ) + + groups = build_periodic_groups( + points, + np.asarray([0.0, 0.0]), + np.asarray([1.0, 1.0]), + [True, False], + tolerance=1.0e-6, + ) + + self.assertEqual(groups, []) + def test_density_string_translation_supports_c_style_blocks(self): density = _compile_density_function( """ @@ -190,11 +356,124 @@ def test_density_string_translation_supports_c_style_blocks(self): self.assertAlmostEqual(density(np.asarray([0.0, 6.0])), 4.0) self.assertAlmostEqual(density(np.asarray([0.0, 0.0])), 5.0) + def test_invalid_density_string_raises_instead_of_falling_back(self): + with self.assertRaises(ValueError): + _compile_density_function("density = ;") + + def test_engine_honors_configured_max_steps_without_python_cap(self): + engine = RelaxationEngine( + self._box_geometry(), + {"parameters": {"controller_step_limit_max": 250, "nr_probes_for_determining_volume": 200}}, + 0.5, + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + [False, False], + ) + + self.assertEqual(engine.max_steps, 250) + + def test_boundary_drift_tolerance_does_not_scale_with_a0(self): + engine = RelaxationEngine( + self._box_geometry(), + {"parameters": {"nr_probes_for_determining_volume": 200}}, + 10.0, + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + [False, False], + ) + + self.assertEqual(engine.boundary_drift_tolerance, BOUNDARY_FUZZ) + + def test_point_change_schedule_uses_legacy_square_steps(self): + engine = RelaxationEngine( + self._box_geometry(), + {"parameters": {"controller_step_limit_max": 100, "nr_probes_for_determining_volume": 200}}, + 0.5, + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + [False, False], + ) + + scheduled_steps = [] + for step in range(1, 30): + engine.step = step + if engine._should_attempt_point_change(): + scheduled_steps.append(step) + + self.assertEqual(scheduled_steps, [11, 14, 19, 26]) + + def test_simply_points_disable_generated_seed_points(self): + simply_points = np.asarray( + [ + [0.1, 0.1], + [0.9, 0.1], + [0.5, 0.8], + ], + dtype=float, + ) + + points, states = _prepare_initial_points( + self._box_geometry(), + 0.25, + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + simply_points, + [False, False], + np.random.default_rng(97), + ) + + self.assertEqual(points.tolist(), simply_points.tolist()) + self.assertEqual(states.tolist(), [STATE_SIMPLE, STATE_SIMPLE, STATE_SIMPLE]) + + def test_user_mobile_points_disable_generated_seed_points(self): + mobile_points = np.asarray( + [ + [0.2, 0.2], + [0.8, 0.2], + [0.5, 0.7], + ], + dtype=float, + ) + + points, _states = _prepare_initial_points( + self._box_geometry(), + 0.25, + np.empty((0, 2), dtype=float), + mobile_points, + np.empty((0, 2), dtype=float), + [False, False], + np.random.default_rng(97), + ) + + self.assertEqual(points.tolist(), mobile_points.tolist()) + + def test_step_limit_waits_for_post_change_settling(self): + engine = RelaxationEngine( + self._box_geometry(), + {"parameters": {"controller_step_limit_max": 20, "nr_probes_for_determining_volume": 200}}, + 0.5, + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + [False, False], + ) + + engine.last_addition_deletion_step = 11 + engine.step = 20 + self.assertFalse(engine._step_limit_reached()) + engine.step = 61 + self.assertTrue(engine._step_limit_reached()) + def test_meshing_is_deterministic_for_identical_inputs(self): kwargs = dict( bounding_box=[[0.0, 0.0], [1.0, 1.0]], objects=[Box([0.1, 0.1], [0.9, 0.9])], a0=0.5, + max_steps=20, + nr_probes_for_determining_volume=500, ) mesh_a = Mesh(**kwargs) diff --git a/tests/nmesh/test_parity.py b/tests/nmesh/test_parity.py new file mode 100644 index 0000000..6efe756 --- /dev/null +++ b/tests/nmesh/test_parity.py @@ -0,0 +1,154 @@ +import pytest + +from nmesh.backend import RawMesh +from nmesh.mesher.parity import ( + assert_canonical_mesh_equal, + canonical_mesh_signature, + compare_mesh_metrics, + mesh_metric_summary, + read_ascii_nmesh, +) + + +def _reference_triangle() -> RawMesh: + return RawMesh( + points=[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]], + simplices=[[0, 1, 2]], + regions=[7], + point_regions=[[7], [7], [7]], + surfaces=[[0, 1], [0, 2], [1, 2]], + links=[(0, 1), (0, 2), (1, 2)], + periodic_point_indices=[[0, 1]], + region_volumes=[0.5], + dim=2, + ) + + +def test_canonical_mesh_signature_compares_topology_after_point_reordering(): + actual = RawMesh( + points=[[1.0, 0.0], [0.0, 1.0], [0.0, 0.0]], + simplices=[[2, 0, 1]], + regions=[7], + point_regions=[[7], [7], [7]], + surfaces=[[2, 0], [1, 2], [0, 1]], + links=[(2, 0), (1, 2), (0, 1)], + periodic_point_indices=[[2, 0]], + region_volumes=[0.5], + dim=2, + ) + + assert_canonical_mesh_equal(actual, _reference_triangle()) + + +def test_canonical_mesh_signature_is_exact_by_default(): + actual = _reference_triangle() + expected = _reference_triangle() + expected.points[1][0] = 1.0 + 1.0e-12 + + with pytest.raises(AssertionError): + assert_canonical_mesh_equal(actual, expected) + + +def test_canonical_mesh_signature_requires_explicit_coordinate_tolerance(): + actual = _reference_triangle() + expected = _reference_triangle() + expected.points[1][0] = 1.0 + 1.0e-12 + + assert_canonical_mesh_equal( + actual, + expected, + coordinate_tolerance=1.0e-10, + ) + + +def test_canonical_mesh_signature_rejects_ambiguous_tolerance(): + mesh = _reference_triangle() + mesh.points[2] = [0.0, 0.5e-9] + + with pytest.raises(ValueError): + canonical_mesh_signature(mesh, coordinate_tolerance=1.0e-9) + + +def test_read_ascii_nmesh_accepts_legacy_surface_and_periodic_rows(tmp_path): + mesh_path = tmp_path / "legacy.nmesh" + mesh_path.write_text( + "\n".join( + [ + "# PYFEM mesh file version 1.0", + "# dim = 2 nodes = 3 simplices = 1 surfaces = 3 periodic = 1", + "3", + "0.0 0.0", + "1.0 0.0", + "0.0 1.0", + "1", + "7 0 1 2", + "3", + "7 -1 0 1", + "7 -1 0 2", + "7 -1 1 2", + "1", + "0 0 1", + "", + ] + ), + encoding="utf-8", + ) + + raw_mesh = read_ascii_nmesh(mesh_path) + + assert raw_mesh.points == [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]] + assert raw_mesh.simplices == [[0, 1, 2]] + assert raw_mesh.regions == [7] + assert raw_mesh.surfaces == [[0, 1], [0, 2], [1, 2]] + assert raw_mesh.periodic_point_indices == [[0, 1]] + assert raw_mesh.links == [(0, 1), (0, 2), (1, 2)] + assert raw_mesh.region_volumes == [0.5] + + +def test_mesh_metric_summary_reports_geometry_quality_metrics(): + summary = mesh_metric_summary(_reference_triangle()) + + assert summary.dim == 2 + assert summary.point_count == 3 + assert summary.simplex_count == 1 + assert summary.region_counts == ((7, 1),) + assert summary.region_volumes == ((7, 0.5),) + assert summary.bbox_min == (0.0, 0.0) + assert summary.bbox_max == (1.0, 1.0) + assert summary.simplex_measure_min == 0.5 + assert summary.simplex_measure_mean == 0.5 + assert summary.simplex_measure_max == 0.5 + + +def test_compare_mesh_metrics_allows_configured_metric_tolerances(): + actual = _reference_triangle() + expected = _reference_triangle() + expected.points[1][0] = 1.02 + expected.region_volumes = [0.51] + + comparison = compare_mesh_metrics( + actual, + expected, + volume_relative_tolerance=0.05, + length_relative_tolerance=0.05, + ) + + assert comparison.passed + + +def test_compare_mesh_metrics_reports_failures_outside_tolerance(): + actual = _reference_triangle() + expected = _reference_triangle() + expected.points[1][0] = 2.0 + expected.region_volumes = [1.0] + + comparison = compare_mesh_metrics( + actual, + expected, + volume_relative_tolerance=0.05, + length_relative_tolerance=0.05, + ) + + assert not comparison.passed + assert any("bbox_max" in failure for failure in comparison.failures) + assert any("region_volumes" in failure for failure in comparison.failures) diff --git a/tests/nmesh_test.py b/tests/nmesh_test.py index 5b0858f..4de5eaa 100644 --- a/tests/nmesh_test.py +++ b/tests/nmesh_test.py @@ -52,7 +52,13 @@ def test_mesh_generation_produces_topology(self): bb = [[0,0,0], [1,1,1]] obj = nmesh.Box([0.2,0.2,0.2], [0.8,0.8,0.8]) - m = nmesh.Mesh(bounding_box=bb, objects=[obj], a0=0.1) + m = nmesh.Mesh( + bounding_box=bb, + objects=[obj], + a0=0.3, + max_steps=20, + nr_probes_for_determining_volume=500, + ) self.assertGreater(len(m.points), 0) self.assertGreater(len(m.simplices), 0) self.assertTrue(all(region == 1 for region in m.regions)) @@ -68,6 +74,8 @@ def test_periodic_mesh_generation(self): a0=0.5, periodic=[True, False], mesh_bounding_box=True, + max_steps=20, + nr_probes_for_determining_volume=500, ) self.assertGreater(len(m.points), 0) @@ -89,6 +97,8 @@ def test_hint_mesh_seeds_generation(self): mesh_bounding_box=True, hints=[(seed, nmesh.Box([0.1, 0.1], [0.9, 0.9]))], a0=0.5, + max_steps=20, + nr_probes_for_determining_volume=500, ) self.assertIn([0.2, 0.2], mesh.points) @@ -107,6 +117,8 @@ def callback(piece, step, mesh_info): objects=[nmesh.Box([0.0, 0.0], [1.0, 1.0])], a0=0.5, callback=(callback, 2), + max_steps=20, + nr_probes_for_determining_volume=500, ) self.assertGreater(len(mesh.points), 0) diff --git a/tools/nmesh_parity_compare.py b/tools/nmesh_parity_compare.py new file mode 100644 index 0000000..2eff79f --- /dev/null +++ b/tools/nmesh_parity_compare.py @@ -0,0 +1,247 @@ +#!/usr/bin/env python3 +"""Run modern nmesh scenarios and compare them with legacy mesh artifacts. + +Examples: + python tools/nmesh_parity_compare.py --write-new --output-dir parity/new + python tools/nmesh_parity_compare.py --legacy-dir parity/legacy + python tools/nmesh_parity_compare.py --legacy-command "legacy_runner {scenario} {output}" +""" + +from __future__ import annotations + +import argparse +import json +import shlex +import subprocess +import tempfile +from collections.abc import Callable +from dataclasses import asdict, dataclass +from pathlib import Path + +import nmesh +from nmesh.backend import RawMesh +from nmesh.mesher.parity import ( + assert_canonical_mesh_equal, + compare_mesh_metrics, + mesh_metric_summary, + read_ascii_nmesh, +) + + +@dataclass(frozen=True, slots=True) +class Scenario: + """One modern-vs-legacy nmesh comparison case.""" + + name: str + build: Callable[[], RawMesh] + + +def _box2d() -> RawMesh: + return nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[nmesh.Box([0.1, 0.1], [0.9, 0.9])], + a0=0.35, + max_steps=20, + nr_probes_for_determining_volume=500, + ).raw_mesh + + +def _box3d() -> RawMesh: + return nmesh.Mesh( + bounding_box=[[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]], + objects=[nmesh.Box([0.2, 0.2, 0.2], [0.8, 0.8, 0.8])], + a0=0.3, + max_steps=20, + nr_probes_for_determining_volume=500, + ).raw_mesh + + +def _adjacent2d() -> RawMesh: + return nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[ + nmesh.Box([0.0, 0.0], [0.5, 1.0]), + nmesh.Box([0.5, 0.0], [1.0, 1.0]), + ], + a0=0.5, + max_steps=20, + nr_probes_for_determining_volume=500, + ).raw_mesh + + +def _concave2d() -> RawMesh: + return nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[ + nmesh.difference( + nmesh.Box([0.0, 0.0], [1.0, 1.0]), + [nmesh.Box([0.35, 0.35], [0.65, 0.65])], + ) + ], + a0=0.25, + max_steps=20, + nr_probes_for_determining_volume=1000, + ).raw_mesh + + +def _periodic2d() -> RawMesh: + return nmesh.Mesh( + bounding_box=[[0.0, 0.0], [1.0, 1.0]], + objects=[nmesh.Box([0.0, 0.0], [1.0, 1.0])], + mesh_bounding_box=True, + periodic=[True, True], + a0=0.4, + max_steps=20, + nr_probes_for_determining_volume=600, + ).raw_mesh + + +SCENARIOS: dict[str, Scenario] = { + scenario.name: scenario + for scenario in ( + Scenario("box2d", _box2d), + Scenario("box3d", _box3d), + Scenario("adjacent2d", _adjacent2d), + Scenario("concave2d", _concave2d), + Scenario("periodic2d", _periodic2d), + ) +} + + +def main() -> int: + """Run configured scenarios and print a JSON parity report.""" + + args = _parse_args() + scenario_names = args.scenario or sorted(SCENARIOS) + output_dir = Path(args.output_dir) if args.output_dir else None + if output_dir is not None: + output_dir.mkdir(parents=True, exist_ok=True) + + with tempfile.TemporaryDirectory(prefix="nmesh-parity-") as temp_dir: + temp_path = Path(temp_dir) + report = [] + exit_code = 0 + for scenario_name in scenario_names: + scenario = SCENARIOS[scenario_name] + modern_mesh = scenario.build() + modern_path = (output_dir or temp_path) / f"{scenario_name}.modern.nmesh" + nmesh.write_mesh(modern_mesh, modern_path) + + scenario_report = { + "scenario": scenario_name, + "modern_path": str(modern_path), + "modern_metrics": asdict(mesh_metric_summary(modern_mesh)), + } + legacy_path = _legacy_path_for(args, scenario_name, temp_path) + if legacy_path is None: + scenario_report["status"] = "modern-written" + report.append(scenario_report) + continue + + legacy_mesh = read_ascii_nmesh(legacy_path) + scenario_report["legacy_path"] = str(legacy_path) + scenario_report["legacy_metrics"] = asdict(mesh_metric_summary(legacy_mesh)) + exact_failure = _exact_failure(modern_mesh, legacy_mesh, args.coordinate_tolerance) + comparison = compare_mesh_metrics( + modern_mesh, + legacy_mesh, + count_relative_tolerance=args.count_relative_tolerance, + volume_relative_tolerance=args.volume_relative_tolerance, + length_relative_tolerance=args.length_relative_tolerance, + absolute_tolerance=args.absolute_tolerance, + ) + scenario_report["exact_match"] = exact_failure is None + if exact_failure is not None: + scenario_report["exact_failure"] = exact_failure + scenario_report["metric_passed"] = comparison.passed + scenario_report["metric_failures"] = list(comparison.failures) + scenario_report["status"] = "passed" if comparison.passed else "failed" + if not comparison.passed: + exit_code = 1 + report.append(scenario_report) + + print(json.dumps(report, indent=2, sort_keys=True)) + return exit_code + + +def _parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--scenario", + action="append", + choices=sorted(SCENARIOS), + help="Scenario to run. May be provided multiple times. Defaults to all scenarios.", + ) + parser.add_argument( + "--legacy-dir", + type=Path, + help="Directory containing legacy '.nmesh' or '.legacy.nmesh' files.", + ) + parser.add_argument( + "--legacy-command", + help=( + "Command template used to generate a legacy mesh. The template may " + "contain {scenario} and {output}; it is split with shlex and run without a shell." + ), + ) + parser.add_argument( + "--output-dir", + type=Path, + help="Directory where modern scenario meshes should be written.", + ) + parser.add_argument("--coordinate-tolerance", type=float) + parser.add_argument("--count-relative-tolerance", type=float, default=0.25) + parser.add_argument("--volume-relative-tolerance", type=float, default=0.05) + parser.add_argument("--length-relative-tolerance", type=float, default=0.10) + parser.add_argument("--absolute-tolerance", type=float, default=1.0e-9) + return parser.parse_args() + + +def _legacy_path_for( + args: argparse.Namespace, + scenario_name: str, + temp_path: Path, +) -> Path | None: + if args.legacy_command: + output_path = temp_path / f"{scenario_name}.legacy.nmesh" + command = args.legacy_command.format( + scenario=scenario_name, + output=str(output_path), + ) + subprocess.run(shlex.split(command), check=True) + return output_path + + if args.legacy_dir is None: + return None + + candidates = [ + args.legacy_dir / f"{scenario_name}.legacy.nmesh", + args.legacy_dir / f"{scenario_name}.nmesh", + ] + for candidate in candidates: + if candidate.exists(): + return candidate + raise FileNotFoundError( + f"No legacy mesh found for {scenario_name!r}; expected one of: " + + ", ".join(str(candidate) for candidate in candidates) + ) + + +def _exact_failure( + actual: RawMesh, + expected: RawMesh, + coordinate_tolerance: float | None, +) -> str | None: + try: + assert_canonical_mesh_equal( + actual, + expected, + coordinate_tolerance=coordinate_tolerance, + ) + except AssertionError as exc: + return str(exc) + return None + + +if __name__ == "__main__": + raise SystemExit(main())