Triangulation

The top-level entry point for building a triangular mesh on a 2D domain.

triangulate

admesh.triangulate

triangulate(domain: Domain | str | 'os.PathLike[str]', *, h_max: float | None = None, h_min: float | None = None, size_field: Callable[[ndarray], ndarray] | None = None, user_contribs: tuple[Callable[[ndarray], ndarray], ...] = (), combine: Callable[[list[ndarray]], ndarray] = np.minimum.reduce, background: str = 'uniform', seed: int | None = None, max_iter: int | None = None, initial_points: 'np.ndarray | None' = None, quality_gate: tuple[float, float] = (0.3, 0.6), ttol: float | None = None, dptol: float | None = None, medial_method: str | None = None) -> Mesh

Generate a triangular mesh on domain.

Parameters:

Name	Type	Description	Default
domain	Domain, str, or os.PathLike	Domain object, file path (TOML/JSON/fort.14), or mesh_id string (from ADMESH-Domains registry if installed).	required
h_max	float or None	Target maximum edge length. If None, defaults to bbox_diagonal / 20.	None
h_min	float or None	Minimum edge length for size field composition.	None
size_field	callable or None	Pre-composed size field function.	None
user_contribs	tuple of callables	User-defined size field contributions.	()
combine	callable	Function to combine multiple size fields (default: np.minimum.reduce).	reduce
background	str	Background grid strategy: 'uniform' (default) or 'octree' (adaptive octree-backed size field from spec-029).	'uniform'
seed	int or None	Random seed for reproducibility.	None
max_iter	int or None	Maximum iterations for mesh generation.	None
quality_gate	tuple[float, float]	Advisory (min_q, mean_q) smoke thresholds. Default: (0.30, 0.60) — an MVP port-sanity floor, NOT a binding quality invariant (#140). Quality is driven by h_min/h_max/g; pass (0.0, 0.0) to disable the gate when knobs lower min quality.	(0.3, 0.6)
ttol	float or None	Relative displacement threshold for Delaunay rebuild. Default: 0.27.	None
dptol	float or None	Interior node movement tolerance for convergence. Default: 2e-3.	None
medial_method	(None, 'grid', 'octree', 'vdt')	Optional channel-width size contribution, added to the other contributions and combined through combine. None (default) adds nothing and leaves the output unchanged. The size is clip(LFS / 2, h_min, h_max) where LFS is the local feature size, the distance to the boundary plus the distance to the medial axis. The methods differ in how the medial axis is found: 'grid' uses the grid method of the 2012 port, 'octree' uses the octree leaf graph, and 'vdt' uses the vector distance transform of Kang and Kubatko (2024, https://doi.org/10.5194/gmd-17-1603-2024) with corner pruning. All three share one grid spacing, h_max / 4 (enlarged to keep the grid under four million cells), and use h_min (default h_max / 100) as the lower size bound. When a method is chosen, the initial lattice spacing of the generator is lowered to the smallest size the contribution takes inside the domain (not below h_min), because the generator thins a lattice of that spacing.	None

Returns:

Type	Description
Mesh	Triangulated mesh with quality metrics and boundaries.

Raises:

Type	Description
ValueError	If domain source cannot be resolved, quality gates fail, or medial_method is not one of the values listed above.
ImportError	If registry lookup is attempted without valence-domains installed.
Adapts the v1 :class:`Domain` onto the faithful-port driver
:func:`admesh.routine.triangulate` without modifying it (Constitution
Principle I). Returns a :class:`Mesh` with per-element quality
populated and boundaries derived from the triangulation.

Source code in src/admesh/api.py

def triangulate(
    domain: Domain | str | "os.PathLike[str]",
    *,
    h_max: float | None = None,
    h_min: float | None = None,
    size_field: Callable[[np.ndarray], np.ndarray] | None = None,
    user_contribs: tuple[Callable[[np.ndarray], np.ndarray], ...] = (),
    combine: Callable[[list[np.ndarray]], np.ndarray] = np.minimum.reduce,
    background: str = "uniform",
    seed: int | None = None,
    max_iter: int | None = None,
    initial_points: "np.ndarray | None" = None,
    # Advisory default, NOT a binding invariant (#140):
    # (0.30, 0.60) is an MVP port-sanity smoke floor, caller-overridable. Mesh
    # quality is hyperparameter-driven (h_min/h_max/g); aggressive ratios
    # legitimately lower min quality. Pass (0.0, 0.0) to disable the gate.
    quality_gate: tuple[float, float] = (0.30, 0.60),
    ttol: float | None = None,
    dptol: float | None = None,
    medial_method: str | None = None,
) -> Mesh:
    """Generate a triangular mesh on ``domain``.

    Parameters
    ----------
    domain : Domain, str, or os.PathLike
        Domain object, file path (TOML/JSON/fort.14), or mesh_id string
        (from ADMESH-Domains registry if installed).
    h_max : float or None
        Target maximum edge length. If None, defaults to bbox_diagonal / 20.
    h_min : float or None
        Minimum edge length for size field composition.
    size_field : callable or None
        Pre-composed size field function.
    user_contribs : tuple of callables
        User-defined size field contributions.
    combine : callable
        Function to combine multiple size fields (default: np.minimum.reduce).
    background : str
        Background grid strategy: 'uniform' (default) or 'octree' (adaptive
        octree-backed size field from spec-029).
    seed : int or None
        Random seed for reproducibility.
    max_iter : int or None
        Maximum iterations for mesh generation.
    quality_gate : tuple[float, float]
        Advisory (min_q, mean_q) smoke thresholds. Default: (0.30, 0.60) —
        an MVP port-sanity floor, NOT a binding quality invariant
        (#140). Quality is driven by h_min/h_max/g;
        pass (0.0, 0.0) to disable the gate when knobs lower min quality.
    ttol : float or None
        Relative displacement threshold for Delaunay rebuild. Default: 0.27.
    dptol : float or None
        Interior node movement tolerance for convergence. Default: 2e-3.
    medial_method : {None, 'grid', 'octree', 'vdt'}
        Optional channel-width size contribution, added to the other
        contributions and combined through ``combine``. ``None`` (default)
        adds nothing and leaves the output unchanged. The size is
        ``clip(LFS / 2, h_min, h_max)`` where ``LFS`` is the local feature
        size, the distance to the boundary plus the distance to the medial
        axis. The methods differ in how the medial axis is found:
        ``'grid'`` uses the grid method of the 2012 port, ``'octree'`` uses the
        octree leaf graph, and ``'vdt'`` uses the vector distance transform
        of Kang and Kubatko (2024, https://doi.org/10.5194/gmd-17-1603-2024)
        with corner pruning. All three share one grid spacing, ``h_max / 4``
        (enlarged to keep the grid under four million cells), and use
        ``h_min`` (default ``h_max / 100``) as the lower size bound. When a
        method is chosen, the initial lattice spacing of the generator is
        lowered to the smallest size the contribution takes inside the domain
        (not below ``h_min``), because the generator thins a lattice of that
        spacing.

    Returns
    -------
    Mesh
        Triangulated mesh with quality metrics and boundaries.

    Raises
    ------
    ValueError
        If domain source cannot be resolved, quality gates fail, or
        ``medial_method`` is not one of the values listed above.
    ImportError
        If registry lookup is attempted without valence-domains installed.

    Adapts the v1 :class:`Domain` onto the faithful-port driver
    :func:`admesh.routine.triangulate` without modifying it (Constitution
    Principle I). Returns a :class:`Mesh` with per-element quality
    populated and boundaries derived from the triangulation.
    """
    if medial_method is not None and medial_method not in _MEDIAL_METHODS:
        raise ValueError(
            f"triangulate: medial_method must be None or one of "
            f"{_MEDIAL_METHODS}, got {medial_method!r}"
        )

    # Lazy imports — keeps `import admesh` cheap and avoids a hard
    # dependency cycle with the faithful-port modules at import time.
    from admesh._stages.domains import Domain as _PortDomain
    from admesh._stages.quality import mesh_quality
    from admesh._stages.routine import triangulate as _routine_triangulate

    # Load domain from file or registry if it's a string; adapt if it's a port Domain.
    api_domain = None  # Track whether we have an api.Domain (may have bc_segments)
    if isinstance(domain, _PortDomain):
        # Input is already a faithful-port Domain; use it directly.
        port_domain = domain
        bbox = domain.bbox
        h0_default = max(_bbox_diag(bbox) / 20.0, 1e-6)
        h0 = float(h_max) if h_max is not None else h0_default
        pfix = np.asarray(domain.fixed_points, dtype=np.float64) if domain.fixed_points is not None else np.empty((0, 2), dtype=np.float64)
    else:
        # Input should be an api.Domain or a file/registry path
        if not isinstance(domain, Domain):
            domain = _load_domain_from_source(domain)

        api_domain = domain  # Track the api.Domain for bc_segments access later
        # Resolve h0 from h_max, falling back to a fraction of the bbox diagonal.
        if h_max is None:
            h0 = max(_bbox_diag(domain.bbox) / 20.0, 1e-6)
        else:
            h0 = float(h_max)

        # Adapter: the faithful-port `Domain` carries the SDF + fixed points.
        pfix = domain.pfix
        if pfix is None:
            pfix = np.empty((0, 2), dtype=np.float64)
        pfix = np.asarray(pfix, dtype=np.float64)
        port_domain = _PortDomain(
            name="api_v1",
            fd=domain.sdf,
            bbox=domain.bbox,
            fixed_points=pfix,
        )

    # Build kwargs for the driver. seed / niter only forwarded when the
    # caller supplied them — let the routine apply its own defaults.
    opts: dict[str, object] = {}
    if seed is not None:
        opts["seed"] = int(seed)
    if max_iter is not None:
        opts["niter"] = int(max_iter)
    if ttol is not None:
        opts["ttol"] = float(ttol)
    if dptol is not None:
        opts["dptol"] = float(dptol)
    if initial_points is not None:
        opts["initial_points"] = np.asarray(initial_points, dtype=np.float64)

    # Resolve the size field. Cases:
    #
    #   1. Caller passed a pre-composed `size_field=`. They've already
    #      done their own composition; we use it as-is. If they ALSO
    #      passed `user_contribs=` we warn — those would be ignored
    #      otherwise, which silently violates the contract.
    #   2. Caller passed `user_contribs=`. Wrap them via
    #      `compose_size_field` with `size_field` (if any) as the sole
    #      Phase-1 builtin. Default combiner is `np.minimum.reduce`.
    #   3. Caller passed h_min or h_max but no user_contribs — auto-create
    #      a uniform clamping size field so the bounds are not silently ignored.
    #   4. Neither — uniform sizing falls through (`fh=None`).
    medial_fn = None
    if medial_method is not None:
        _m_hmax = float(h0)
        _m_hmin = float(h_min) if h_min is not None else _m_hmax / 100.0
        medial_fn, _m_floor = _medial_contribution(
            medial_method,
            getattr(port_domain, "fd"),
            tuple(float(b) for b in port_domain.bbox),
            h_max=_m_hmax,
            h_min=_m_hmin,
        )
        # The driver treats ``h0`` as the smallest target edge length (it seeds
        # a lattice of that spacing and thins it by the size field), so a
        # contribution that refines below ``h_max`` needs a finer ``h0``.
        h0 = min(h0, max(_m_floor, _m_hmin))

    if size_field is not None and user_contribs:
        warnings.warn(
            "triangulate: both `size_field` and `user_contribs` were "
            "supplied; ignoring `user_contribs` (the pre-composed "
            "`size_field` already encodes its own composition).",
            UserWarning,
            stacklevel=2,
        )
        if medial_fn is None:
            fh = size_field
        else:
            from admesh.size_field import compose_size_field

            fh = compose_size_field(
                builtins=(size_field,),
                user_contribs=(medial_fn,),
                combine=combine,
                hmin=h_min,
                hmax=h_max,
            )
    elif user_contribs or h_min is not None or h_max is not None or medial_fn is not None:
        from admesh.size_field import compose_size_field

        # Build Phase-1 builtins: include size_field if provided
        builtins_phase1 = (size_field,) if size_field is not None else ()

        # Build Phase-2 user contributions
        user_phase2 = tuple(user_contribs)
        _has_other = bool(builtins_phase1) or bool(user_phase2)

        # If h_min/h_max specified without other size fields, create a bounded field.
        # NOTE (#65 / spec 025 Step 3 DEFERRED): wiring build_h() here as the
        # unconditional default degrades MVP convex-domain min_q 0.30 -> 0.22
        # (curvature+medial size-gradient -> low-quality distmesh transition
        # tris), violating the constitutional MVP quality gate + spec 025
        # AC-005/AC-006. Deferred pending operator decision (make conditional on
        # domain features / tune scales / revise gate). Steps 1+2 (Domain.bathymetry
        # + from_mesh extraction) shipped — additive, no behavior change.
        if not _has_other:
            # Create a clamping size field directly to avoid warning noise.
            def _clamped_uniform(pts):
                pts_arr = np.asarray(pts, dtype=np.float64)
                n = pts_arr.shape[0]
                # Return uniform field, already clamped to [h_min, h_max]
                result = np.full(n, h_max if h_max is not None else h_min or 1.0, dtype=np.float64)
                if h_min is not None and h_max is not None:
                    result[:] = h_max  # Use h_max as the target (conservative for refinement)
                return result
            user_phase2 = (_clamped_uniform,)

        if medial_fn is not None:
            user_phase2 = (*user_phase2, medial_fn)

        fh = compose_size_field(
            builtins=builtins_phase1,
            user_contribs=user_phase2,
            combine=combine,
            hmin=h_min,
            hmax=h_max,
        )
    elif size_field is None and (h_min is not None or h_max is not None):
        # Fix #37: h_min/h_max were silently ignored when no user_contribs
        # provided. Auto-compose a uniform field that clamps to [h_min, h_max].
        from admesh.size_field import compose_size_field

        _h_max_val = h_max if h_max is not None else h0
        _uniform = lambda pts: np.full(len(pts), _h_max_val, dtype=float)  # noqa: E731
        fh = compose_size_field(
            builtins=(_uniform,),
            user_contribs=(),
            combine=combine,
            hmin=h_min,
            hmax=h_max,
        )
    else:
        fh = size_field  # may still be None — uniform sizing

    if background == "octree":
        from admesh.octree import octree_size_field
        _bbox = port_domain.bbox
        _hmax_o = float(h_max) if h_max is not None else h0
        _hmin_o = float(h_min) if h_min is not None else _hmax_o / 100.0
        _base = fh if fh is not None else (
            lambda pts: np.full(len(np.atleast_2d(pts)), _hmax_o, dtype=float))
        class _BBoxShim:
            bbox = _bbox
        fh = octree_size_field(_BBoxShim, _base, h_min=_hmin_o, h_max=_hmax_o)
    elif background != "uniform":
        raise ValueError(f"triangulate: background must be 'uniform' or 'octree', got {background!r}")

    # Issue #2: if the domain carries explicit boundary vertices (pts), seed
    # intermediate points along each edge so short boundary segments get
    # adequate coverage even when the 2-D lattice is coarse.
    # Use getattr because admesh.domains.Domain (MVP class) lacks `pts`.
    _pts = getattr(domain, "pts", None)
    domain_pts = _pts if _pts is not None else getattr(domain, "boundary_polygon", None)
    if domain_pts is not None:
        boundary_seeds = _seed_boundary_1d(
            np.asarray(domain_pts, dtype=np.float64), fh, h0
        )
        if boundary_seeds.size:
            pfix = (
                np.vstack([pfix, boundary_seeds]) if pfix.size else boundary_seeds
            )
            # _PortDomain uses `fd`; api.Domain uses `sdf` — handle both.
            _sdf = getattr(domain, "sdf", None) or getattr(domain, "fd", None)
            port_domain = _PortDomain(
                name="api_v1",
                fd=_sdf,
                bbox=domain.bbox,
                fixed_points=pfix,
            )

    p, t = _routine_triangulate(port_domain, h0=h0, fh=fh, **opts)
    nodes = np.asarray(p, dtype=np.float64)
    elements = np.asarray(t, dtype=np.int64)

    min_q, mean_q, q_per = mesh_quality(nodes, elements)
    gate_min, gate_mean = quality_gate
    if min_q < gate_min:
        raise ValueError(
            f"triangulate: min_q {min_q:.3f} < quality_gate[0] {gate_min:.2f}"
        )
    if mean_q < gate_mean:
        raise ValueError(
            f"triangulate: mean_q {mean_q:.3f} < quality_gate[1] {gate_mean:.2f}"
        )

    # Derive boundary segments. If the caller pre-declared bc_segments
    # on the Domain, pass them through verbatim — they're the user's
    # contract about the output. Otherwise default-label every closed
    # boundary ring as MAINLAND.
    if api_domain and api_domain.bc_segments:
        # Validate that all node_ids in bc_segments are valid for the new mesh.
        # When a Domain comes from Domain.from_mesh(old_mesh), bc_segments may
        # reference node ids from the old mesh, which are invalid for the new
        # triangulation (e.g., old mesh had 9933 boundary nodes, new has 103).
        max_node_id = -1
        valid = True
        for seg in api_domain.bc_segments:
            if seg.node_ids.size > 0:
                seg_max = int(np.max(seg.node_ids))
                max_node_id = max(max_node_id, seg_max)
                if seg_max >= len(nodes):
                    valid = False
                    break

        if valid and max_node_id >= 0:
            boundaries = tuple(api_domain.bc_segments)
        else:
            # Fall back to deriving boundaries from the triangulation
            if max_node_id >= len(nodes):
                warnings.warn(
                    f"triangulate: Domain.bc_segments reference node ids outside the generated mesh "
                    f"(max id {max_node_id}, n_nodes {len(nodes)}); deriving boundaries from the triangulation instead.",
                    UserWarning,
                    stacklevel=2
                )
            boundaries = _derive_boundary_segments(elements, nodes)
    else:
        boundaries = _derive_boundary_segments(elements, nodes)

    return Mesh(
        nodes=nodes,
        elements=elements,
        boundaries=boundaries,
        bathymetry=None,
        quality=np.asarray(q_per, dtype=np.float64),
        title="",
    )

triangulate_batch

Run triangulate over many domains on a process pool. Results come back in input order and match a serial loop.

admesh.triangulate_batch

triangulate_batch(domains: Sequence, *, n_jobs: int | None = None, **kwargs) -> list[Mesh]

Triangulate multiple domains on a process pool.

Parameters:

Name	Type	Description	Default
domains	sequence	Sequence of Domain objects, file paths, or registry slugs — anything :func:triangulate accepts.	required
n_jobs	int or None	Maximum number of worker processes. If None, defaults to min(len(domains), os.cpu_count() or 1). If 1, runs sequentially in-process without spawning workers. If < 1, raises ValueError.	None
**kwargs		Keyword arguments forwarded unchanged to every :func:triangulate call (h_max, h_min, size_field, seed, max_iter, quality_gate, etc.).	{}

Returns:

Type	Description
list of Mesh	Meshes in input order. Results from sequential in-process execution (n_jobs=1) are numerically identical to pool results; the pool is an execution detail only.

Raises:

Type	Description
ValueError	If n_jobs is not None and < 1.
TypeError	If any domain or kwargs are not picklable (raised before pool creation to avoid silent serialization failures in workers). Message includes the domain index (for domains) or 'kwargs'.

Notes

When n_jobs=1 or len(domains) == 1, the function runs a sequential loop without spawning a process pool. This avoids pickling overhead and permits unpicklable domains such as those with lambda SDF functions.

For parallel execution (n_jobs > 1), each domain is pickled before the pool is created. Domains with lambda or local-scope SDF functions will raise TypeError. Use registry slugs, file paths, or module-level SDF callables instead, or set n_jobs=1.

Source code in src/admesh/batch.py

def triangulate_batch(
    domains: Sequence,
    *,
    n_jobs: int | None = None,
    **kwargs,
) -> list[Mesh]:
    """Triangulate multiple domains on a process pool.

    Parameters
    ----------
    domains : sequence
        Sequence of Domain objects, file paths, or registry slugs — anything
        :func:`triangulate` accepts.
    n_jobs : int or None, optional
        Maximum number of worker processes. If None, defaults to
        ``min(len(domains), os.cpu_count() or 1)``.
        If 1, runs sequentially in-process without spawning workers.
        If < 1, raises ValueError.
    **kwargs
        Keyword arguments forwarded unchanged to every :func:`triangulate`
        call (h_max, h_min, size_field, seed, max_iter, quality_gate, etc.).

    Returns
    -------
    list of Mesh
        Meshes in input order. Results from sequential in-process
        execution (n_jobs=1) are numerically identical to pool results;
        the pool is an execution detail only.

    Raises
    ------
    ValueError
        If n_jobs is not None and < 1.
    TypeError
        If any domain or kwargs are not picklable (raised before pool
        creation to avoid silent serialization failures in workers).
        Message includes the domain index (for domains) or 'kwargs'.

    Notes
    -----
    When n_jobs=1 or ``len(domains) == 1``, the function runs a
    sequential loop without spawning a process pool. This avoids
    pickling overhead and permits unpicklable domains such as those
    with lambda SDF functions.

    For parallel execution (n_jobs > 1), each domain is pickled before
    the pool is created. Domains with lambda or local-scope SDF
    functions will raise TypeError. Use registry slugs, file paths,
    or module-level SDF callables instead, or set n_jobs=1.
    """
    domains_list = list(domains)

    # Resolve n_jobs
    if n_jobs is None:
        n_jobs = min(len(domains_list), os.cpu_count() or 1)
    if n_jobs < 1:
        raise ValueError(f"n_jobs must be >= 1, got {n_jobs}")

    # Empty input
    if not domains_list:
        return []

    # Sequential execution for n_jobs=1 or single domain
    if n_jobs == 1 or len(domains_list) == 1:
        return [triangulate(domain, **kwargs) for domain in domains_list]

    # Parallel execution: pickle-check all domains before pool creation
    for i, domain in enumerate(domains_list):
        try:
            pickle.dumps(domain)
        except (TypeError, AttributeError, pickle.PicklingError) as e:
            raise TypeError(
                f"Domain at index {i} is not picklable and cannot be sent to "
                f"worker processes. Suggestions:\n"
                f"  - Use a module-level SDF callable instead of a lambda\n"
                f"  - Use a file path (.toml, .json, .14) or registry slug\n"
                f"  - Set n_jobs=1 to run sequentially in-process"
            ) from e

    # Pickle-check kwargs once
    try:
        pickle.dumps(kwargs)
    except (TypeError, AttributeError, pickle.PicklingError) as e:
        raise TypeError(
            "kwargs are not picklable and cannot be sent to worker processes. "
            "Suggestions:\n"
            "  - Avoid callables with lambda or local scope\n"
            "  - Use module-level callables for size_field, bathymetry, etc.\n"
            "  - Set n_jobs=1 to run sequentially in-process"
        ) from e

    # Create pool and submit tasks
    mp_context = multiprocessing.get_context("spawn")
    with ProcessPoolExecutor(max_workers=n_jobs, mp_context=mp_context) as executor:
        futures = [
            executor.submit(_worker_triangulate, domain, kwargs)
            for domain in domains_list
        ]
        results = [future.result() for future in futures]

    return results