Validation¶

def magnet_source(label: str, center_m: Vec3, size_m: Vec3, domain: DomainConfig, *, at: Mapping[str, Quantity] | None = None) -> MagnetSource | None:
    """Build the `MagnetSource` for one body, or None if it is not a magnet; a magnet with no usable `Magnetization` direction raises `ValueError`.

    Args:
        label: The body's label, carrying its tags.
        center_m: The body's center, in metres.
        size_m: The body's axis-aligned extent, in metres.
        domain: The domain config whose materials and tags the label matches.
        at: The operating point the remanence is evaluated at; None is room temperature.
    """
    material: MaterialProperties | None = material_for_label(label, domain.materials)
    if material is None or not material.is_magnet:
        return None

    point: Mapping[str, Quantity] = at if at is not None else _room_temperature()
    magnitude: float | None = material.magnetic.magnetization_magnitude(at=point)  # |M| in A/m

    assert magnitude is not None  # is_magnet (remanence present) guarantees a magnitude; narrow the Optional

    magnetization: Magnetization | None = magnetization_for_label(label, domain.tags)

    if magnetization is None or magnetization.direction is None:
        raise ValueError(f"magnet body {label!r} (material {material.tag!r}) has no usable Magnetization direction tag")

    polarization: Vec3 = magnetization.direction * (magnitude * PERMEABILITY_OF_VACUUM)  # J = μ₀·|M| (tesla)

    return MagnetSource(center=center_m, size=size_m, polarization=polarization)

def sources_from_assembly(
    assembly: Compound, domain: DomainConfig, *, at: Mapping[str, Quantity] | None = None
) -> list[MagnetSource]:  # pragma: no cover - reads a live build123d assembly (OCC, no stubs); proven by the gated integration test
    """Walk a built build123d assembly and return a `MagnetSource` per magnet leaf, reading only its label and bounding box (mm → m).

    Args:
        assembly: The built assembly, whose magnets turn only by multiples of 90° and so stay axis-aligned.
        domain: The domain config whose materials and tags the leaf labels match.
        at: The operating point the remanence is evaluated at; None is room temperature.
    """
    sources: list[MagnetSource] = []

    for leaf in assembly.leaves:  # pyright: ignore[reportUnknownMemberType, reportUnknownVariableType]  # build123d ships no stubs
        leaf: Part = leaf  # pyright: ignore[ reportUnknownVariableType]  # build123d ships no stubs
        box = leaf.located(leaf.global_location).bounding_box()  # pyright: ignore[reportUnknownMemberType, reportUnknownVariableType, reportUnknownArgumentType]
        center = box.center()  # pyright: ignore[reportUnknownMemberType, reportUnknownVariableType]
        size = box.size  # pyright: ignore[reportUnknownMemberType, reportUnknownVariableType]
        center_m: Vec3 = Vec3(float(center.X), float(center.Y), float(center.Z)) * MM_TO_M  # pyright: ignore[reportUnknownArgumentType, reportUnknownMemberType]
        size_m: Vec3 = Vec3(float(size.X), float(size.Y), float(size.Z)) * MM_TO_M  # pyright: ignore[reportUnknownArgumentType, reportUnknownMemberType]

        source: MagnetSource | None = magnet_source(str(leaf.label), center_m, size_m, domain, at=at)  # pyright: ignore[reportUnknownArgumentType, reportUnknownMemberType]
        if source is not None:
            sources.append(source)

    return sources

def material_for_label(label: str, materials: list[MaterialProperties]) -> MaterialProperties | None:
    """The material whose `tag` equals (casefold) the label's `Material` tag value, or None — mirrors `MeshEntity._match_material`."""
    value: str | None = parse_tags(label).get(TagKey.MATERIAL)
    if value is None:
        return None
    return next((material for material in materials if material.tag.casefold() == value.casefold()), None)

def magnetization_for_label(label: str, tags: list[EntityTag]) -> Magnetization | None:
    """The usable (directed) `Magnetization` carried by an `EntityTag` matching the label's tags, or None — mirrors `elmer.sim`'s magnet wiring."""
    matched: list[EntityTag] = [t for t in tags if t.matches(parse_tags(label))]
    conditions = conditions_for(matched, Physics.MAGNETOSTATICS, ConditionTarget.BODY)
    magnetizations: list[Magnetization] = [c for c in conditions if isinstance(c, Magnetization)]

    return next((m for m in magnetizations if m.direction is not None), None)

def gap_line(sources: list[MagnetSource], samples: int) -> list[Vec3]:
    """
    A line of `samples` points along Z through the centroid of the magnet centres — for a dual-Halbach array it threads the air gap.

    Returns a single mid-point when `samples < 2`; raises `ValueError` for an empty list.
    """
    if not sources:
        raise ValueError("gap_line needs at least one source")

    n: int = len(sources)
    cx: float = sum(s.center.x for s in sources) / n
    cy: float = sum(s.center.y for s in sources) / n
    zmin: float = min(s.center.z - s.size.z / 2 for s in sources)
    zmax: float = max(s.center.z + s.size.z / 2 for s in sources)

    if samples < 2:
        return [Vec3(cx, cy, (zmin + zmax) / 2)]

    step: float = (zmax - zmin) / (samples - 1)
    return [Vec3(cx, cy, zmin + step * i) for i in range(samples)]

def sample_magpylib_b(
    sources: list[MagnetSource], points: np.typing.NDArray[np.float64]
) -> np.typing.NDArray[np.float64]:  # pragma: no cover - needs magpylib (viz extra); proven by the gated integration test
    """Evaluate the analytical **B** (tesla) of `sources` at each point: one SI magpylib `Cuboid` per source, superposed via `Collection.getB`."""
    import magpylib as magpy
    import numpy as np
    from magpylib.magnet import Cuboid

    cuboids: list[Cuboid] = [
        Cuboid(
            position=(s.center.x, s.center.y, s.center.z),
            dimension=(s.size.x, s.size.y, s.size.z),
            polarization=(s.polarization.x, s.polarization.y, s.polarization.z),
        )
        for s in sources
    ]
    # getB superposes every source's field at every observer point; errstate silences the rare on-edge singularity
    # (those points are dropped as NaN by the caller before comparison).
    with np.errstate(invalid="ignore", divide="ignore"):
        b = magpy.Collection(*cuboids).getB(points)
    # getB squeezes a single observer down to shape (3,); normalize to (N, 3) so callers always get rows.
    return cast("np.typing.NDArray[np.float64]", np.atleast_2d(b))

def sample_vtu_b(vtu_path: str | Path, points: list[Vec3]) -> list[Vec3]:  # pragma: no cover - needs pyvista (viz extra) + a solved case.vtu; proven by the gated integration test
    """Interpolate the solved Elmer **B** (tesla) at `points` from a `case.vtu` (pyvista `sample`; elemental arrays are projected to point data first)."""
    import numpy as np
    import pyvista as pv

    grid: pv.UnstructuredGrid = cast("pv.UnstructuredGrid", pv.read(str(vtu_path)))
    b_name: str = select_flux_density_array(list(grid.array_names))
    sampled: pv.DataSet = cast("pv.DataSet", grid.cell_data_to_point_data() if b_name in grid.cell_data else grid)

    point_array = np.array([[p.x, p.y, p.z] for p in points], dtype=float)
    probe: pv.PolyData = pv.PolyData(point_array)
    result: pv.PolyData = cast("pv.PolyData", probe.sample(sampled))

    b_values = np.atleast_2d(np.asarray(result.point_data[b_name], dtype=float))
    rows: list[list[float]] = cast("list[list[float]]", b_values.tolist())
    return [Vec3(row[0], row[1], row[2]) for row in rows]

def extract_vtu_b(
    vtu_path: str | Path,
) -> tuple[np.typing.NDArray[np.float64], np.typing.NDArray[np.float64]]:  # pragma: no cover - needs pyvista (viz extra) + a solved case.vtu; proven by the gated integration test
    """Extract the solved Elmer **B** (tesla) from a `case.vtu` as `(positions, B)` arrays of shape `(n_points, 3)`."""
    import numpy as np
    import pyvista as pv

    grid: pv.UnstructuredGrid = cast("pv.UnstructuredGrid", pv.read(str(vtu_path)))

    b_name: str = select_flux_density_array(list(grid.array_names))

    sampled: pv.DataSet = cast("pv.DataSet", grid.cell_data_to_point_data() if b_name in grid.cell_data else grid)

    positions = np.asarray(sampled.points, dtype=np.float64)
    b_field = np.asarray(sampled.point_data[b_name], dtype=np.float64)

    return positions, b_field

def compare_fields(points: NDArray[np.float64], b_reference: NDArray[np.float64], b_test: NDArray[np.float64]) -> FieldComparison:
    """Compare two `(N, 3)` fields sampled at the same `points` and return a `FieldComparison`; raises `ValueError` on mismatched or empty inputs.

    Args:
        points: The `(N, 3)` sample points, in metres.
        b_reference: The reference field at `points`; points where it is exactly zero are left out of the relative errors.
        b_test: The field under test at `points`.
    """
    if not (len(points) == len(b_reference) == len(b_test)):
        raise ValueError(f"compare_fields needs equal-length inputs, got points={len(points)} reference={len(b_reference)} test={len(b_test)}")
    if len(points) == 0:
        raise ValueError("compare_fields needs at least one sample point")

    n: int = len(points)
    abs_error: NDArray[np.float64] = _row_norms(b_test - b_reference)
    ref_norm: NDArray[np.float64] = _row_norms(b_reference)

    max_index: int = int(np.argmax(abs_error))
    # Cast the row to a typed array first (numpy scalar indexing reads as Any), then `.tolist()` gives python floats.
    worst_row: NDArray[np.float64] = cast("NDArray[np.float64]", points[max_index])
    worst_x, worst_y, worst_z = cast("list[float]", worst_row.tolist())
    max_error_point: Vec3 = Vec3(worst_x, worst_y, worst_z)
    max_abs_error: float = float(cast("float", abs_error[max_index]))

    rms_error: float = _rms(abs_error)
    ref_rms: float = _rms(ref_norm)
    test_rms: float = _rms(_row_norms(b_test))

    # Relative metrics over points with a non-zero reference (a near-null |B_ref| has no meaningful relative error).
    positive: NDArray[np.bool_] = ref_norm > 0
    mean_relative_error: float | None = None
    median_relative_error: float | None = None
    if bool(positive.any()):
        relative: NDArray[np.float64] = cast("NDArray[np.float64]", abs_error[positive] / ref_norm[positive])
        mean_relative_error = float(cast("float", np.mean(relative)))
        median_relative_error = float(cast("float", np.median(relative)))
    rel_rms_error: float | None = rms_error / ref_rms if ref_rms > 0 else None

    return FieldComparison(
        n_points=n,
        rms_error=rms_error,
        max_abs_error=max_abs_error,
        max_error_point=max_error_point,
        ref_rms=ref_rms,
        test_rms=test_rms,
        rel_rms_error=rel_rms_error,
        median_relative_error=median_relative_error,
        mean_relative_error=mean_relative_error,
    )

def grid_spec(
    bounds: Bounds,
    samples: int,
    mode: GridBoundsMode = GridBoundsMode.OUTER,
) -> tuple[Point3, Point3, tuple[int, int, int]]:
    """The `(origin, spacing, dimensions)` of a `pv.ImageData` sampling `bounds` with near-cubic cells; raises `ValueError` for fewer than 8 samples or a zero-volume box.

    Args:
        bounds: The target box, as ``(xmin, xmax, ymin, ymax, zmin, zmax)``.
        samples: The approximate total point count.
        mode: Keep the sampled extent inside the bounds (`INNER`) or covering them (`OUTER`).
    """
    if samples < 8:
        raise ValueError("grid_spec needs at least 8 samples (2 per axis × 3 axes)")
    xmin, xmax, ymin, ymax, zmin, zmax = bounds
    ex: float = xmax - xmin
    ey: float = ymax - ymin
    ez: float = zmax - zmin
    if ex <= 0 or ey <= 0 or ez <= 0:
        raise ValueError("grid_spec requires a non-degenerate bounding box")

    spacing: float = math.pow(ex * ey * ez / samples, 1.0 / 3.0)

    quantize = math.floor if mode == GridBoundsMode.INNER else math.ceil

    nx: int = max(2, quantize(ex / spacing))
    ny: int = max(2, quantize(ey / spacing))
    nz: int = max(2, quantize(ez / spacing))

    return (xmin, ymin, zmin), (spacing, spacing, spacing), (nx, ny, nz)

def bounds_of_sources(sources: list[MagnetSource], *, padding: Vec3) -> Bounds:
    """The axis-aligned bounding box (metres) enclosing every source's full extent, grown by `padding` on each side; raises `ValueError` for no sources."""
    if not sources:
        raise ValueError("bounds_of_sources needs at least one source")
    xmin: float = min(s.center.x - s.size.x / 2 for s in sources)
    xmax: float = max(s.center.x + s.size.x / 2 for s in sources)
    ymin: float = min(s.center.y - s.size.y / 2 for s in sources)
    ymax: float = max(s.center.y + s.size.y / 2 for s in sources)
    zmin: float = min(s.center.z - s.size.z / 2 for s in sources)
    zmax: float = max(s.center.z + s.size.z / 2 for s in sources)
    return (xmin - padding.x, xmax + padding.x, ymin - padding.y, ymax + padding.y, zmin - padding.z, zmax + padding.z)

def bounded_grid(bounds: Bounds, samples: int) -> pv.ImageData:  # pragma: no cover - needs pyvista (viz extra); proven by the gated integration test
    """A `pv.ImageData` grid spanning `bounds` with ~cubic cells and `samples` in total, ready to carry a sampled vector field as point data."""
    import pyvista as pv

    origin, spacing, dims = grid_spec(bounds, samples)
    return pv.ImageData(dimensions=dims, spacing=spacing, origin=origin)

def magpylib_field_grid(
    sources: list[MagnetSource], *, padding: Vec3 = DEFAULT_PADDING, samples: int = 32
) -> pv.ImageData:  # pragma: no cover - needs magpylib + pyvista (viz extra); proven by the gated integration test
    """Sample the analytical **B** of `sources` onto a regular grid over their padded bounding box, a `pv.ImageData` ready for the streamline path.

    Args:
        sources: The magnets.
        padding: The air added around the magnets on each side, in metres.
        samples: The approximate total point count.
    """
    import numpy as np
    import pyvista as pv

    bounds: Bounds = bounds_of_sources(sources, padding=padding)
    origin, spacing, dims = grid_spec(bounds, samples)
    grid: pv.ImageData = pv.ImageData(dimensions=dims, spacing=spacing, origin=origin)
    # pyvista's `.points` has an Unknown dtype; cast pins it so the magpylib sampler gets a typed (N, 3) array.
    points: np.typing.NDArray[np.float64] = cast("np.typing.NDArray[np.float64]", grid.points)
    grid.point_data[FLUX_DENSITY] = sample_magpylib_b(sources, points)
    return grid

def load_gmsh_mesh(
    msh_path: str | Path,
) -> pv.UnstructuredGrid:  # pragma: no cover - needs pyvista + meshio (viz extra); proven by the gated integration test
    """Load a Gmsh `.msh` (via meshio) into a `pv.UnstructuredGrid` with the region names attached as the `"Region_Name"` cell-data array."""
    import meshio
    import numpy as np
    import pyvista as pv

    msh_data = cast("meshio.Mesh", meshio.read(str(msh_path)))  # pyright: ignore[reportUnknownMemberType]  # meshio.read returns a Mesh, but meshio ships no stubs

    cells = cast("list[meshio.CellBlock]", msh_data.cells)
    total_cells = sum(len(cast("np.ndarray", cell_block.data)) for cell_block in cells)

    cell_names = np.empty(total_cells, dtype=object)

    # Map the per-block cell_sets masks back to global cell index positions.
    for set_name, cell_set_masks in msh_data.cell_sets.items():
        if set_name == "gmsh:bounding_entities":  # gmsh metadata set, not a region
            continue

        global_idx_offset = 0
        for block_idx, mask in enumerate(cast("list[np.ndarray | None]", cell_set_masks)):
            block_length = len(cast("np.ndarray", cells[block_idx].data))

            if mask is not None and len(mask) > 0:
                # Mask can be a boolean array or an array of local integer indices.
                if mask.dtype == bool:
                    local_indices = np.where(mask)[0]
                else:
                    local_indices = mask

                global_indices = local_indices + global_idx_offset
                cell_names[global_indices] = set_name

            global_idx_offset += block_length

    # Clear cell_sets before conversion — the PyVista reader crashes on it.
    msh_data.cell_sets = {}
    mesh = pv.from_meshio(msh_data)

    mesh.cell_data["Region_Name"] = cell_names

    return mesh