"""Container holding about the geometry and medium properties common to all types of simulations."""
from __future__ import annotations
from typing import TYPE_CHECKING, Any
import autograd.numpy as np
from pydantic import Field, field_validator
if TYPE_CHECKING:
from typing import Literal
from pydantic import NonNegativeInt
from tidy3d.compat import Self
from tidy3d.components.material.types import StructureMediumType
from .data.utils import CustomSpatialDataType
from .grid.grid import Grid
from .types import (
Ax,
Bound,
Coordinate,
InterpMethod,
PermittivityComponent,
PlotScale,
Shapely,
Size,
)
from .viz import PlotParams
from tidy3d.components.material.tcad.charge import (
ChargeConductorMedium,
SemiconductorMedium,
)
from tidy3d.components.material.tcad.heat import (
AnisotropicConductivity,
SolidMedium,
SolidSpec,
)
from tidy3d.components.material.types import MultiPhysicsMediumType3D
from tidy3d.components.tcad.doping import (
ConstantDoping,
CustomDoping,
DopingBoxType,
GaussianDoping,
)
from tidy3d.components.tcad.doping_utils import aggregate_doping_seam_aware
from tidy3d.components.tcad.viz import HEAT_SOURCE_CMAP
from tidy3d.constants import CONDUCTIVITY, THERMAL_CONDUCTIVITY, inf
from tidy3d.exceptions import SetupError, Tidy3dError
from tidy3d.log import log
from .base import Tidy3dBaseModel, cached_property
from .data.utils import (
SpatialDataArray,
TetrahedralGridDataset,
TriangularGridDataset,
UnstructuredGridDataset,
)
from .geometry.base import Box
from .geometry.utils import merging_geometries_on_plane
from .grid.grid import Coords
from .material.multi_physics import MultiPhysicsMedium
from .medium import (
AbstractCustomMedium,
AbstractMedium,
AbstractPerturbationMedium,
Medium,
Medium2D,
)
from .structure import Structure
from .types import TYPE_TAG_STR, LengthUnit, PriorityMode
from .validators import assert_unique_names
from .viz import (
MEDIUM_CMAP,
STRUCTURE_EPS_CMAP,
STRUCTURE_EPS_CMAP_R,
STRUCTURE_HEAT_COND_CMAP,
add_ax_if_none,
equal_aspect,
plot_params_fluid,
plot_params_structure,
plot_scene_3d,
polygon_path,
)
# Sentinel to indicate matplotlib hasn't been loaded yet
_MPL_NOT_LOADED = object()
# Module-level mpl reference for lazy loading and test mocking
# - _MPL_NOT_LOADED: not yet imported (will import on first use)
# - None: explicitly unavailable (for testing fallback behavior)
# - module: matplotlib module when successfully imported
mpl: Any = _MPL_NOT_LOADED
# maximum number of mediums supported
MAX_NUM_MEDIUMS = 65530
# # maximum geometry count in a single structure
# MAX_GEOMETRY_COUNT = 5000
# warn and error out if the same medium is present in too many structures
WARN_STRUCTURES_PER_MEDIUM = 200
MAX_STRUCTURES_PER_MEDIUM = 1_000
def _get_colormap(reverse: bool = False) -> str:
return STRUCTURE_EPS_CMAP_R if reverse else STRUCTURE_EPS_CMAP
def _heat_conductivity_scalar(conductivity: float | AnisotropicConductivity) -> float:
"""Scalar for grayscale plotting: mean of principals for a tensor, pass-through for a scalar."""
if isinstance(conductivity, AnisotropicConductivity):
return (conductivity.xx + conductivity.yy + conductivity.zz) / 3
return conductivity
[docs]
class Scene(Tidy3dBaseModel):
"""Contains generic information about the geometry and medium properties common to all types of
simulations.
Example
-------
>>> sim = Scene(
... structures=[
... Structure(
... geometry=Box(size=(1, 1, 1), center=(0, 0, 0)),
... medium=Medium(permittivity=2.0),
... ),
... ],
... medium=Medium(permittivity=3.0),
... )
"""
medium: MultiPhysicsMediumType3D = Field(
default_factory=Medium,
title="Background Medium",
description="Background medium of scene, defaults to vacuum if not specified.",
discriminator=TYPE_TAG_STR,
)
structures: tuple[Structure, ...] | None = Field(
(),
title="Structures",
description="Tuple of structures present in scene. "
"Note: In regions of spatial overlap between structures, "
"material properties are dictated by structure of higher priority. "
"The priority for structure of `priority=None` is set automatically "
"based on `structure_priority_mode`. For structures of equal priority, "
"the structure added later to the structure list takes precedence.",
)
structure_priority_mode: PriorityMode = Field(
"equal",
title="Structure Priority Setting",
description="This field only affects structures of `priority=None`. "
"If `equal`, the priority of those structures is set to 0; if `conductor`, "
"the priority of structures made of `LossyMetalMedium` is set to 90, "
"`PECMedium` to 100, and others to 0.",
)
plot_length_units: LengthUnit | None = Field(
"μm",
title="Plot Units",
description="When set to a supported ``LengthUnit``, "
"plots will be produced with proper scaling of axes and "
"include the desired unit specifier in labels.",
)
""" Validating setup """
_unique_structure_names = assert_unique_names("structures")
@field_validator("structures")
@classmethod
def _validate_mediums(cls, val: tuple[Structure, ...] | None) -> tuple[Structure, ...] | None:
"""Error if too many mediums present. Warn if different mediums have the same name."""
if val is None:
return val
mediums = {structure.medium for structure in val}
if len(mediums) > MAX_NUM_MEDIUMS:
raise SetupError(
f"Tidy3D only supports {MAX_NUM_MEDIUMS} distinct mediums."
f"{len(mediums)} were supplied."
)
medium_names = [medium.name for medium in mediums if medium.name is not None]
if len(medium_names) != len(set(medium_names)):
log.warning(
"Different mediums with the same name were detected. "
"This may error in future Tidy3D versions, and using unique names for distinct "
"media is recommended."
)
return val
# @field_validator("structures")
# @classmethod
# def _validate_num_geometries(cls, val):
# """Error if too many geometries in a single structure."""
# if val is None:
# return val
# for i, structure in enumerate(val):
# for geometry in flatten_groups(structure.geometry, flatten_transformed=True):
# count = sum(
# 1
# for g in traverse_geometries(geometry)
# if not isinstance(g, (GeometryGroup, ClipOperation))
# )
# if count > MAX_GEOMETRY_COUNT:
# raise SetupError(
# f"Structure at 'structures[{i}]' has {count} geometries that cannot be "
# f"flattened. A maximum of {MAX_GEOMETRY_COUNT} is supported due to "
# f"preprocessing performance."
# )
# return val
@field_validator("structures")
@classmethod
def _validate_structures_per_medium(
cls, val: tuple[Structure, ...] | None
) -> tuple[Structure, ...] | None:
"""Error if too many structures share the same medium; suggest using GeometryGroup."""
if val is None:
return val
# if total structures are <= warn limit, the constraint cannot be violated.
if len(val) <= WARN_STRUCTURES_PER_MEDIUM:
return val
# Count structures per medium
counts = {}
get = counts.get
for structure in val:
key = structure.medium
new_count = get(key, 0) + 1
# Exit early to avoid slow counting if many structures present
if new_count > MAX_STRUCTURES_PER_MEDIUM:
raise SetupError(
f"More than {MAX_STRUCTURES_PER_MEDIUM} structures use the same medium. "
"For performance, use a 'GeometryGroup' or boolean operations to combine "
"geometries that share a medium."
)
counts[key] = new_count
# Now check if we should warn
for count in counts.values():
if count > WARN_STRUCTURES_PER_MEDIUM:
log.warning(
f"More than {WARN_STRUCTURES_PER_MEDIUM} structures use the same medium. "
"For performance, use a 'GeometryGroup' or boolean operations to combine "
"geometries that share a medium."
)
break
return val
""" Accounting """
@cached_property
def bounds(self) -> Bound:
"""Automatically defined scene's bounds based on present structures. Infinite dimensions
are ignored. If the scene contains no structures, the bounds are set to
(-1, -1, -1), (1, 1, 1). Similarly, if along a given axis all structures extend infinitely,
the bounds along that axis are set from -1 to 1.
Returns
-------
tuple[float, float, float], tuple[float, float, float]
Min and max bounds packaged as ``(minx, miny, minz), (maxx, maxy, maxz)``.
"""
bounds = tuple(structure.geometry.bounds for structure in self.structures)
return (
tuple(min((b[i] for b, _ in bounds if b[i] != -inf), default=-1) for i in range(3)),
tuple(max((b[i] for _, b in bounds if b[i] != inf), default=1) for i in range(3)),
)
@cached_property
def size(self) -> Size:
"""Automatically defined scene's size.
Returns
-------
tuple[float, float, float]
Scene's size.
"""
return tuple(bmax - bmin for bmin, bmax in zip(self.bounds[0], self.bounds[1]))
@cached_property
def center(self) -> Coordinate:
"""Automatically defined scene's center.
Returns
-------
tuple[float, float, float]
Scene's center.
"""
return tuple(0.5 * (bmin + bmax) for bmin, bmax in zip(self.bounds[0], self.bounds[1]))
@cached_property
def box(self) -> Box:
"""Automatically defined scene's :class:`.Box`.
Returns
-------
Box
Scene's box.
"""
return Box(center=self.center, size=self.size)
@cached_property
def mediums(self) -> set[StructureMediumType]:
"""Returns set of distinct :class:`.AbstractMedium` in scene.
Returns
-------
list[:class:`.AbstractMedium`]
Set of distinct mediums in the scene.
"""
medium_dict = {self.medium: None}
medium_dict.update({structure.medium: None for structure in self.structures})
return list(medium_dict.keys())
@cached_property
def medium_map(self) -> dict[StructureMediumType, NonNegativeInt]:
"""Returns dict mapping medium to index in material.
``medium_map[medium]`` returns unique global index of :class:`.AbstractMedium` in scene.
Returns
-------
dict[:class:`.AbstractMedium`, int]
Mapping between distinct mediums to index in scene.
"""
return {medium: index for index, medium in enumerate(self.mediums)}
@cached_property
def sorted_structures(self) -> list[Structure]:
"""Returns a list of sorted structures based on their priority.In the sorted list,
latter added structures take higher priority.
Returns
-------
List[:class:`.Structure`]
"""
return Structure._sort_structures(self.structures, self.structure_priority_mode)
@cached_property
def background_structure(self) -> Structure:
"""Returns structure representing the background of the :class:`.Scene`."""
geometry = Box(size=(inf, inf, inf))
return Structure(geometry=geometry, medium=self.medium)
@cached_property
def all_structures(self) -> list[Structure]:
"""List of all structures in the simulation including the background."""
return [self.background_structure, *self.sorted_structures]
[docs]
@staticmethod
def intersecting_structures(
test_object: Box, structures: tuple[Structure, ...]
) -> tuple[Structure, ...]:
"""From a given list of structures, returns a list of :class:`.Structure` that intersect
with the ``test_object``, if it is a surface, or its surfaces, if it is a volume.
Parameters
-------
test_object : :class:`.Box`
Object for which intersecting media are to be detected.
structures : list[:class:`.AbstractMedium`]
List of structures whose media will be tested.
Returns
-------
list[:class:`.Structure`]
Set of distinct structures that intersect with the given surface, or with the surfaces
of the given volume.
"""
if test_object.size.count(0.0) == 1:
# get all merged structures on the test_object, which is already planar
normal_axis_index = test_object.size.index(0.0)
dim = "xyz"[normal_axis_index]
pos = test_object.center[normal_axis_index]
xyz_kwargs = {dim: pos}
structures_merged = []
for structure in structures:
intersections = structure.geometry.intersections_plane(**xyz_kwargs)
if len(intersections) > 0:
structures_merged.append(structure)
return structures_merged
# if the test object is a volume, test each surface recursively
surfaces = test_object.surfaces_with_exclusion(**test_object.model_dump())
structures_merged = []
for surface in surfaces:
structures_merged += Scene.intersecting_structures(surface, structures)
return structures_merged
""" Plotting General """
@staticmethod
def _get_plot_lims(
bounds: Bound,
x: float | None = None,
y: float | None = None,
z: float | None = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
) -> tuple[tuple[float, float], tuple[float, float]]:
# if no hlim and/or vlim given, the bounds will then be the usual pml bounds
axis, _ = Box.parse_xyz_kwargs(x=x, y=y, z=z)
_, (hmin, vmin) = Box.pop_axis(bounds[0], axis=axis)
_, (hmax, vmax) = Box.pop_axis(bounds[1], axis=axis)
# account for unordered limits
if hlim is None:
hlim = (hmin, hmax)
if vlim is None:
vlim = (vmin, vmax)
if hlim[0] > hlim[1]:
raise Tidy3dError("Error: 'hmin' > 'hmax'")
if vlim[0] > vlim[1]:
raise Tidy3dError("Error: 'vmin' > 'vmax'")
if hlim[0] == hlim[1]:
margin = 0.1 * abs(hlim[0]) if hlim[0] != 0 else 0.05
hlim = (hlim[0] - margin, hlim[1] + margin)
if vlim[0] == vlim[1]:
margin = 0.1 * abs(vlim[0]) if vlim[0] != 0 else 0.05
vlim = (vlim[0] - margin, vlim[1] + margin)
return hlim, vlim
[docs]
@equal_aspect
@add_ax_if_none
def plot(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
ax: Ax = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
fill_structures: bool = True,
**patch_kwargs: Any,
) -> Ax:
"""Plot each of scene's components on a plane defined by one nonzero x,y,z coordinate.
Parameters
----------
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
ax : matplotlib.axes._subplots.Axes = None
Matplotlib axes to plot on, if not specified, one is created.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
fill_structures : bool = True
Whether to fill structures with color or just draw outlines.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
hlim, vlim = Scene._get_plot_lims(bounds=self.bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
ax = self.plot_structures(ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, fill=fill_structures)
ax = self._set_plot_bounds(bounds=self.bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
return ax
[docs]
@equal_aspect
@add_ax_if_none
def plot_structures(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
ax: Ax = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
fill: bool = True,
) -> Ax:
"""Plot each of scene's structures on a plane defined by one nonzero x,y,z coordinate.
Parameters
----------
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
ax : matplotlib.axes._subplots.Axes = None
Matplotlib axes to plot on, if not specified, one is created.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
fill : bool = True
Whether to fill structures with color or just draw outlines.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
medium_shapes = self._get_structures_2dbox(
structures=self.to_static().sorted_structures, x=x, y=y, z=z, hlim=hlim, vlim=vlim
)
medium_map = self.medium_map
for medium, shape in medium_shapes:
mat_index = medium_map[medium]
ax = self._plot_shape_structure(
medium=medium,
mat_index=mat_index,
shape=shape,
ax=ax,
fill=fill,
)
# clean up the axis display
axis, _ = Box.parse_xyz_kwargs(x=x, y=y, z=z)
ax = self.box.add_ax_lims(axis=axis, ax=ax)
ax = self._set_plot_bounds(bounds=self.bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
# Add the default axis labels, tick labels, and title
ax = Box.add_ax_labels_and_title(
ax=ax, x=x, y=y, z=z, plot_length_units=self.plot_length_units
)
return ax
def _plot_shape_structure(
self,
medium: MultiPhysicsMediumType3D,
mat_index: int,
shape: Shapely,
ax: Ax,
fill: bool = True,
) -> Ax:
"""Plot a structure's cross section shape for a given medium."""
plot_params_struct = self._get_structure_plot_params(
medium=medium,
mat_index=mat_index,
fill=fill,
)
ax = self.box.plot_shape(shape=shape, plot_params=plot_params_struct, ax=ax)
return ax
def _get_structure_plot_params(
self,
mat_index: int,
medium: MultiPhysicsMediumType3D,
fill: bool = True,
) -> PlotParams:
"""Constructs the plot parameters for a given medium in scene.plot()."""
plot_params = plot_params_structure.copy(update={"linewidth": 0})
if isinstance(medium, MultiPhysicsMedium):
is_pec = medium.optical is not None and medium.optical.is_pec
is_pmc = medium.optical is not None and medium.optical.is_pmc
is_time_modulated = medium.optical is not None and medium.optical.is_time_modulated
else:
is_pec = medium.is_pec
is_pmc = medium.is_pmc
is_time_modulated = medium.is_time_modulated
if mat_index == 0 or medium == self.medium:
# background medium
plot_params = plot_params.copy(update={"facecolor": "white", "edgecolor": "white"})
elif is_pec:
# perfect electrical conductor
plot_params = plot_params.copy(
update={"facecolor": "gold", "edgecolor": "k", "linewidth": 1}
)
elif is_pmc:
# perfect magnetic conductor
plot_params = plot_params.copy(
update={"facecolor": "purple", "edgecolor": "k", "linewidth": 1}
)
elif is_time_modulated:
# time modulated medium
plot_params = plot_params.copy(
update={"facecolor": "red", "linewidth": 0, "hatch": "x*"}
)
elif isinstance(medium, Medium2D):
# 2d material
plot_params = plot_params.copy(update={"edgecolor": "k", "linewidth": 1})
elif isinstance(medium, Medium):
# regular medium
facecolor = MEDIUM_CMAP[(mat_index - 1) % len(MEDIUM_CMAP)]
plot_params = plot_params.copy(update={"facecolor": facecolor})
if hasattr(medium, "viz_spec"):
if medium.viz_spec is not None:
plot_params = plot_params.override_with_viz_spec(medium.viz_spec)
else:
# regular medium
facecolor = MEDIUM_CMAP[(mat_index - 1) % len(MEDIUM_CMAP)]
plot_params = plot_params.copy(update={"facecolor": facecolor})
if hasattr(medium, "viz_spec"):
if medium.viz_spec is not None:
plot_params = plot_params.override_with_viz_spec(medium.viz_spec)
if not fill:
plot_params = plot_params.copy(update={"fill": False})
if plot_params.linewidth == 0:
plot_params = plot_params.copy(update={"linewidth": 1})
return plot_params
@staticmethod
def _add_cbar(
vmin: float,
vmax: float,
label: str,
cmap: str,
ax: Ax = None,
norm: mpl.colors.Normalize | None = None,
) -> None:
"""Add a colorbar to plot."""
import matplotlib as mpl
import matplotlib.pylab as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.15)
if norm is None:
norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
mappable = mpl.cm.ScalarMappable(norm=norm, cmap=cmap)
plt.colorbar(mappable, cax=cax, label=label)
@staticmethod
def _set_plot_bounds(
bounds: Bound,
ax: Ax,
x: float | None = None,
y: float | None = None,
z: float | None = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
) -> Ax:
"""Sets the xy limits of the scene at a plane, useful after plotting.
Parameters
----------
ax : matplotlib.axes._subplots.Axes
Matplotlib axes to set bounds on.
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
Returns
-------
matplotlib.axes._subplots.Axes
The axes after setting the boundaries.
"""
hlim, vlim = Scene._get_plot_lims(bounds=bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
ax.set_xlim(hlim)
ax.set_ylim(vlim)
return ax
def _get_structures_2dbox(
self,
structures: list[Structure],
x: float | None = None,
y: float | None = None,
z: float | None = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
) -> list[tuple[Medium, Shapely]]:
"""Compute list of shapes to plot on 2d box specified by (x_min, x_max), (y_min, y_max).
Parameters
----------
structures : list[:class:`.Structure`]
list of structures to filter on the plane.
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
Returns
-------
list[tuple[:class:`.AbstractMedium`, shapely.geometry.base.BaseGeometry]]
List of shapes and mediums on the plane.
"""
# if no hlim and/or vlim given, the bounds will then be the usual pml bounds
axis, _ = Box.parse_xyz_kwargs(x=x, y=y, z=z)
_, (hmin, vmin) = Box.pop_axis(self.bounds[0], axis=axis)
_, (hmax, vmax) = Box.pop_axis(self.bounds[1], axis=axis)
if hlim is not None:
(hmin, hmax) = hlim
if vlim is not None:
(vmin, vmax) = vlim
# get center and size with h, v
h_center = (hmin + hmax) / 2.0
v_center = (vmin + vmax) / 2.0
h_size = (hmax - hmin) or inf
v_size = (vmax - vmin) or inf
axis, center_normal = Box.parse_xyz_kwargs(x=x, y=y, z=z)
center = Box.unpop_axis(center_normal, (h_center, v_center), axis=axis)
size = Box.unpop_axis(0.0, (h_size, v_size), axis=axis)
plane = Box(center=center, size=size)
medium_shapes = []
for structure in structures:
intersections = plane.intersections_with(structure.geometry, section_tolerance_2d=True)
for shape in intersections:
if not shape.is_empty:
shape = Box.evaluate_inf_shape(shape)
medium_shapes.append((structure.medium, shape))
return medium_shapes
@staticmethod
def _filter_structures_plane_medium(
structures: list[Structure], plane: Box, section_tolerance_2d: bool = False
) -> list[tuple[Medium, Shapely]]:
"""Compute list of shapes to plot on plane. Overlaps are removed or merged depending on
medium.
Parameters
----------
structures : list[:class:`.Structure`]
List of structures to filter on the plane.
plane : Box
Plane specification.
section_tolerance_2d : bool = False
If ``True``, pass the 2D section tolerance through to the underlying geometry
section queries.
Returns
-------
list[tuple[:class:`.AbstractMedium`, shapely.geometry.base.BaseGeometry]]
List of shapes and mediums on the plane after merging.
"""
medium_list = [structure.medium for structure in structures]
return Scene._filter_structures_plane(
structures=structures,
plane=plane,
property_list=medium_list,
section_tolerance_2d=section_tolerance_2d,
)
@staticmethod
def _filter_structures_plane(
structures: list[Structure],
plane: Box,
property_list: list[Any],
section_tolerance_2d: bool = False,
) -> list[tuple[Medium, Shapely]]:
"""Compute list of shapes to plot on plane. Overlaps are removed or merged depending on
provided property_list.
Parameters
----------
structures : list[:class:`.Structure`]
List of structures to filter on the plane.
plane : Box
Plane specification.
property_list : List = None
Property value for each structure.
section_tolerance_2d : bool = False
If ``True``, pass the 2D section tolerance through to the underlying geometry
section queries.
Returns
-------
list[tuple[:class:`.AbstractMedium`, shapely.geometry.base.BaseGeometry]]
List of shapes and their property value on the plane after merging.
"""
return merging_geometries_on_plane(
[structure.geometry for structure in structures],
plane,
property_list,
section_tolerance_2d=section_tolerance_2d,
)
""" Plotting Optical """
[docs]
@equal_aspect
@add_ax_if_none
def plot_eps(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
freq: float | None = None,
alpha: float | None = None,
ax: Ax = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
eps_lim: tuple[float | None, float | None] = (None, None),
scale: PlotScale = "lin",
) -> Ax:
"""Plot each of scene's components on a plane defined by one nonzero x,y,z coordinate.
The permittivity is plotted in grayscale based on its value at the specified frequency.
Parameters
----------
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
freq : float = None
Frequency to evaluate the relative permittivity of all mediums.
If not specified, evaluates at infinite frequency.
alpha : float = None
Opacity of the structures being plotted.
Defaults to the structure default alpha.
ax : matplotlib.axes._subplots.Axes = None
Matplotlib axes to plot on, if not specified, one is created.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
eps_lim : Tuple[float, float] = None
Custom limits for eps coloring.
scale : PlotScale = "lin"
Scale for the plot. Either 'lin' for linear, 'log' for log10, 'symlog' for symmetric logarithmic (linear near zero, logarithmic elsewhere), or 'dB' for decibel scale.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
hlim, vlim = Scene._get_plot_lims(bounds=self.bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
ax = self.plot_structures_eps(
freq=freq,
cbar=True,
alpha=alpha,
ax=ax,
x=x,
y=y,
z=z,
hlim=hlim,
vlim=vlim,
eps_lim=eps_lim,
scale=scale,
)
ax = self._set_plot_bounds(bounds=self.bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
return ax
[docs]
@equal_aspect
@add_ax_if_none
def plot_structures_eps(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
freq: float | None = None,
alpha: float | None = None,
cbar: bool = True,
reverse: bool = False,
eps_lim: tuple[float | None, float | None] = (None, None),
scale: PlotScale = "lin",
ax: Ax = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
grid: Grid = None,
eps_component: PermittivityComponent | None = None,
) -> Ax:
"""Plot each of scene's structures on a plane defined by one nonzero x,y,z coordinate.
The permittivity is plotted in grayscale based on its value at the specified frequency.
Parameters
----------
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
freq : float = None
Frequency to evaluate the relative permittivity of all mediums.
If not specified, evaluates at infinite frequency.
reverse : bool = False
If ``False``, the highest permittivity is plotted in black.
If ``True``, it is plotteed in white (suitable for black backgrounds).
cbar : bool = True
Whether to plot a colorbar for the relative permittivity.
alpha : float = None
Opacity of the structures being plotted.
Defaults to the structure default alpha.
eps_lim : tuple[float, float] = None
Custom limits for eps coloring.
scale : PlotScale = "lin"
Scale for the plot. Either 'lin' for linear, 'log' for log10, 'symlog' for symmetric logarithmic (linear near zero, logarithmic elsewhere), or 'dB' for decibel scale.
ax : matplotlib.axes._subplots.Axes = None
Matplotlib axes to plot on, if not specified, one is created.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
eps_component : Optional[PermittivityComponent] = None
Component of the permittivity tensor to plot for anisotropic materials,
e.g. ``"xx"``, ``"yy"``, ``"zz"``, ``"xy"``, ``"yz"``, ...
Defaults to ``None``, which returns the average of the diagonal values.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
return self.plot_structures_property(
x=x,
y=y,
z=z,
freq=freq,
alpha=alpha,
cbar=cbar,
reverse=reverse,
limits=eps_lim,
scale=scale,
ax=ax,
hlim=hlim,
vlim=vlim,
grid=grid,
property="eps",
eps_component=eps_component,
)
[docs]
@equal_aspect
@add_ax_if_none
def plot_structures_property(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
freq: float | None = None,
alpha: float | None = None,
cbar: bool = True,
reverse: bool = False,
limits: tuple[float | None, float | None] = (None, None),
scale: PlotScale = "lin",
ax: Ax = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
grid: Grid = None,
property: Literal["eps", "doping", "N_a", "N_d"] = "eps",
eps_component: PermittivityComponent | None = None,
) -> Ax:
"""Plot each of scene's structures on a plane defined by one nonzero x,y,z coordinate.
The permittivity is plotted in grayscale based on its value at the specified frequency.
Parameters
----------
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
freq : float = None
Frequency to evaluate the relative permittivity of all mediums.
If not specified, evaluates at infinite frequency.
reverse : bool = False
If ``False``, the highest permittivity is plotted in black.
If ``True``, it is plotteed in white (suitable for black backgrounds).
cbar : bool = True
Whether to plot a colorbar for the relative permittivity.
alpha : float = None
Opacity of the structures being plotted.
Defaults to the structure default alpha.
limits : tuple[float, float] = None
Custom coloring limits for the property to plot.
scale : PlotScale = "lin"
Scale for the plot. Either 'lin' for linear, 'log' for log10, or 'dB' for decibel scale.
For log scale with negative values, the absolute value is taken before log transformation.
ax : matplotlib.axes._subplots.Axes = None
Matplotlib axes to plot on, if not specified, one is created.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
property: Literal["eps", "doping", "N_a", "N_d"] = "eps"
Indicates the property to plot for the structures. Currently supported properties
are ["eps", "doping", "N_a", "N_d"]
eps_component : Optional[PermittivityComponent] = None
Component of the permittivity tensor to plot for anisotropic materials,
e.g. ``"xx"``, ``"yy"``, ``"zz"``, ``"xy"``, ``"yz"``, ...
Defaults to ``None``, which returns the average of the diagonal values.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
import matplotlib as mpl
structures = self.sorted_structures
# alpha is None just means plot without any transparency
if alpha is None:
alpha = 1
if alpha <= 0:
return ax
need_filtered_shaped = False
if property == "eps":
need_filtered_shaped = alpha < 1 and not isinstance(self.medium, AbstractCustomMedium)
if property in ["N_d", "N_a", "doping"]:
need_filtered_shaped = alpha < 1
if need_filtered_shaped:
axis, position = Box.parse_xyz_kwargs(x=x, y=y, z=z)
center = Box.unpop_axis(position, (0, 0), axis=axis)
size = Box.unpop_axis(0, (inf, inf), axis=axis)
plane = Box(center=center, size=size)
# for doping background structure could be a non-doping structure
# that needs to be rendered
if property in ["N_d", "N_a", "doping"]:
structures = [self.background_structure, *list(structures)]
medium_shapes = self._filter_structures_plane_medium(
structures=structures, plane=plane, section_tolerance_2d=True
)
else:
structures = [self.background_structure, *list(structures)]
medium_shapes = self._get_structures_2dbox(
structures=structures, x=x, y=y, z=z, hlim=hlim, vlim=vlim
)
property_min, property_max = limits
if property_min is None or property_max is None:
if property == "eps":
eps_min_sim, eps_max_sim = self.eps_bounds(freq=freq, eps_component=eps_component)
property_min = property_min if property_min is not None else eps_min_sim
property_max = property_max if property_max is not None else eps_max_sim
linthresh = 1e-2 * property_min
if property in ["N_d", "N_a", "doping"]:
acceptor_limits, donor_limits = self.doping_bounds()
acceptor_abs_min, donor_abs_min = self.doping_absolute_minimum()
if property == "N_d":
property_min = property_min if property_min is not None else donor_limits[0]
property_max = property_max if property_max is not None else donor_limits[1]
linthresh = donor_abs_min
elif property == "N_a":
property_min = property_min if property_min is not None else acceptor_limits[0]
property_max = property_max if property_max is not None else acceptor_limits[1]
linthresh = acceptor_abs_min
elif property == "doping":
property_min = property_min if property_min is not None else -donor_limits[1]
property_max = property_max if property_max is not None else acceptor_limits[1]
linthresh = min(acceptor_abs_min, donor_abs_min)
if np.isclose(linthresh, 0.0):
# fallback to default linthresh of 1e-3
linthresh = 1e-3
else:
if np.isclose(property_min, 0.0) or property_min < 0.0:
linthresh = 1e-3
else:
linthresh = 1e-2 * np.abs(property_min)
if scale == "lin":
norm = mpl.colors.Normalize(vmin=property_min, vmax=property_max)
elif scale == "log":
# LogNorm doesn't work with negative values, so we need to handle this case
if property_min <= 0 or property_max <= 0:
raise SetupError(
f"Log scale cannot be used with non-positive values. "
f"Property range: [{property_min}, {property_max}]. "
f"Consider using 'symlog' scale instead."
)
norm = mpl.colors.LogNorm(vmin=property_min, vmax=property_max)
elif scale == "symlog":
norm = mpl.colors.SymLogNorm(linthresh=linthresh, vmin=property_min, vmax=property_max)
else:
raise SetupError(
f"The scale '{scale}' is not supported for plotting structures property."
)
for medium, shape in medium_shapes:
if property in ["doping", "N_a", "N_d"]:
if not isinstance(medium.charge, SemiconductorMedium):
ax = self._plot_shape_structure_heat_charge_property(
alpha=alpha,
medium=medium,
property_val_min=property_min,
property_val_max=property_max,
reverse=reverse,
shape=shape,
ax=ax,
property="doping",
)
else:
self._pcolormesh_shape_doping_box(
x,
y,
z,
alpha,
medium,
shape,
ax,
property,
norm,
)
else:
# if the background medium is custom medium, it needs to be rendered separately
if medium == self.medium and need_filtered_shaped:
continue
# no need to add patches for custom medium
if not isinstance(medium, AbstractCustomMedium):
ax = self._plot_shape_structure_eps(
freq=freq,
alpha=alpha,
medium=medium,
eps_min=property_min,
eps_max=property_max,
reverse=reverse,
shape=shape,
ax=ax,
eps_component=eps_component,
norm=norm,
)
else:
# For custom medium, apply pcolormesh clipped by the shape.
self._pcolormesh_shape_custom_medium_structure_eps(
x,
y,
z,
freq,
alpha,
medium,
property_min,
property_max,
reverse,
shape,
ax,
grid,
eps_component=eps_component,
norm=norm,
)
if cbar:
if property in ["doping", "N_a", "N_d"]:
Scene._add_cbar(
vmin=property_min,
vmax=property_max,
label=r"$\rm{Doping} \#/cm^3$",
cmap=HEAT_SOURCE_CMAP,
ax=ax,
norm=norm,
)
else:
self._add_cbar_eps(
eps_min=property_min, eps_max=property_max, ax=ax, reverse=reverse, norm=norm
)
# clean up the axis display
axis, _ = Box.parse_xyz_kwargs(x=x, y=y, z=z)
ax = self.box.add_ax_lims(axis=axis, ax=ax)
ax = self._set_plot_bounds(bounds=self.bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
# Add the default axis labels, tick labels, and title
ax = Box.add_ax_labels_and_title(
ax=ax, x=x, y=y, z=z, plot_length_units=self.plot_length_units
)
return ax
@staticmethod
def _add_cbar_eps(
eps_min: float,
eps_max: float,
ax: Ax = None,
reverse: bool = False,
norm: mpl.colors.Normalize | None = None,
) -> None:
"""Add a permittivity colorbar to plot."""
Scene._add_cbar(
vmin=eps_min,
vmax=eps_max,
label=r"$\epsilon_r$",
cmap=_get_colormap(reverse=reverse),
ax=ax,
norm=norm,
)
@staticmethod
def _eps_bounds(
medium_list: list[Medium],
freq: float | None = None,
eps_component: PermittivityComponent | None = None,
) -> tuple[float, float]:
"""Compute range of (real) permittivity present in the mediums at frequency "freq"."""
medium_list = [medium for medium in medium_list if not medium.is_pec]
eps_list = [medium._eps_plot(freq, eps_component) for medium in medium_list]
eps_list = [eps for eps in eps_list if eps is not None]
eps_min = min(eps_list, default=1)
eps_max = max(eps_list, default=1)
# custom medium, the min and max in the supplied dataset over all components and
# spatial locations.
for mat in [medium for medium in medium_list if isinstance(medium, AbstractCustomMedium)]:
mat_epsmin, mat_epsmax = mat._eps_bounds(frequency=freq, eps_component=eps_component)
eps_min = min(eps_min, mat_epsmin)
eps_max = max(eps_max, mat_epsmax)
return eps_min, eps_max
[docs]
def eps_bounds(
self, freq: float | None = None, eps_component: str | None = None
) -> tuple[float, float]:
"""Compute range of (real) permittivity present in the scene at frequency "freq".
Parameters
----------
freq : float = None
Frequency to evaluate the relative permittivity of all mediums.
If not specified, evaluates at infinite frequency.
eps_component : Optional[PermittivityComponent] = None
Component of the permittivity tensor to plot for anisotropic materials,
e.g. ``"xx"``, ``"yy"``, ``"zz"``, ``"xy"``, ``"yz"``, ...
Defaults to ``None``, which returns the average of the diagonal values.
Returns
-------
tuple[float, float]
Minimal and maximal values of relative permittivity in scene.
"""
medium_list = [self.medium, *list(self.mediums)]
return self._eps_bounds(medium_list=medium_list, freq=freq, eps_component=eps_component)
def _pcolormesh_shape_custom_medium_structure_eps(
self,
x: float,
y: float,
z: float,
freq: float,
alpha: float,
medium: Medium,
eps_min: float,
eps_max: float,
reverse: bool,
shape: Shapely,
ax: Ax,
grid: Grid,
eps_component: PermittivityComponent | None = None,
norm: mpl.colors.Normalize | None = None,
) -> None:
"""
Plot shape made of custom medium with ``pcolormesh``.
"""
coords = "xyz"
normal_axis_ind, normal_position = Box.parse_xyz_kwargs(x=x, y=y, z=z)
normal_axis, plane_axes = Box.pop_axis(coords, normal_axis_ind)
comp2ind = {dim + dim: index for dim, index in zip("xyz", range(3))}
# make grid for eps interpolation
# we will do this by combining shape bounds and points where custom eps is provided
shape_bounds = shape.bounds
rmin, rmax = [*shape_bounds[:2]], [*shape_bounds[2:]]
rmin.insert(normal_axis_ind, normal_position)
rmax.insert(normal_axis_ind, normal_position)
if grid is None:
plane_axes_inds = [0, 1, 2]
plane_axes_inds.pop(normal_axis_ind)
eps_diag = medium.eps_dataarray_freq(frequency=freq)
# handle unstructured data case
if isinstance(eps_diag[0], UnstructuredGridDataset):
if (
isinstance(eps_diag[0], TriangularGridDataset)
and eps_diag[0].normal_axis != normal_axis_ind
):
# if we trying to visualize 2d unstructured data not along its normal direction
# we need to extract line slice that lies in the visualization plane
# note that after this eps_diag[] will be SpatialDataArray's
eps_diag = list(eps_diag)
for dim in range(3):
eps_diag[dim] = eps_diag[dim].plane_slice(
axis=normal_axis_ind, pos=normal_position
)
else:
# Select the permittivity component to plot
if eps_component in comp2ind:
eps = eps_diag[comp2ind[eps_component]]
else:
# default to plotting the mean of the diagonal elements
eps = (eps_diag[0] + eps_diag[1] + eps_diag[2]) / 3
if isinstance(eps, TetrahedralGridDataset):
# extract slice if volumetric unstructured data
eps = eps.plane_slice(axis=normal_axis_ind, pos=normal_position)
eps = eps.real
# at this point eps_mean is TriangularGridDataset and we just plot it directly
# with applying shape mask
cmap_name = _get_colormap(reverse=reverse)
eps.plot(
grid=False,
ax=ax,
cbar=False,
cmap=cmap_name,
vmin=eps_min,
vmax=eps_max,
pcolor_kwargs={
"clip_path": (polygon_path(shape), ax.transData),
"clip_box": ax.bbox,
"alpha": alpha,
},
)
return
# in case when different components of custom medium are defined on different grids
# we will combine all points along each dimension
if (
eps_diag[0].coords == eps_diag[1].coords
and eps_diag[0].coords == eps_diag[2].coords
):
coords_to_insert = [eps_diag[0].coords]
else:
coords_to_insert = [eps_diag[0].coords, eps_diag[1].coords, eps_diag[2].coords]
# actual combining of points along each of plane dimensions
plane_coord = []
for ind, comp in zip(plane_axes_inds, plane_axes):
# first start with an array made of shapes bounds
axis_coords = np.array([rmin[ind], rmax[ind]])
# now add points in between them
for coords in coords_to_insert:
comp_axis_coords = coords[comp]
inds_inside_shape = np.where(
np.logical_and(comp_axis_coords > rmin[ind], comp_axis_coords < rmax[ind])
)[0]
if len(inds_inside_shape) > 0:
axis_coords = np.concatenate(
(axis_coords, comp_axis_coords[inds_inside_shape])
)
# remove duplicates
axis_coords = np.unique(axis_coords)
plane_coord.append(axis_coords)
else:
span_inds = grid.discretize_inds(Box.from_bounds(rmin=rmin, rmax=rmax), extend=True)
# filter negative or too large inds
n_grid = [len(grid_comp) for grid_comp in grid.boundaries.to_list]
span_inds = [
(max(fmin, 0), min(fmax, n_grid[f_ind]))
for f_ind, (fmin, fmax) in enumerate(span_inds)
]
# assemble the coordinate in the 2d plane
plane_coord = []
for plane_axis in range(2):
ind_axis = "xyz".index(plane_axes[plane_axis])
plane_coord.append(grid.boundaries.to_list[ind_axis][slice(*span_inds[ind_axis])])
# prepare `Coords` for interpolation
coord_dict = {
plane_axes[0]: plane_coord[0],
plane_axes[1]: plane_coord[1],
normal_axis: [normal_position],
}
coord_shape = Coords(**coord_dict)
# interpolate permittivity and pick the component to plot
eps_shape = medium.eps_diagonal_on_grid(frequency=freq, coords=coord_shape)
if eps_component in comp2ind:
eps_shape = eps_shape[comp2ind[eps_component]]
else:
eps_shape = np.mean(eps_shape, axis=0)
# remove the normal_axis and take real part
eps_shape = eps_shape.real.mean(axis=normal_axis_ind)
# pcolormesh
plane_xp, plane_yp = np.meshgrid(plane_coord[0], plane_coord[1], indexing="ij")
cmap_name = _get_colormap(reverse=reverse)
ax.pcolormesh(
plane_xp,
plane_yp,
eps_shape,
clip_path=(polygon_path(shape), ax.transData),
cmap=cmap_name,
alpha=alpha,
clip_box=ax.bbox,
norm=norm,
)
@staticmethod
def _get_structure_eps_plot_params(
medium: Medium,
freq: float,
eps_min: float,
eps_max: float,
reverse: bool = False,
alpha: float | None = None,
eps_component: PermittivityComponent | None = None,
norm: mpl.colors.Normalize | None = None,
) -> PlotParams:
"""Constructs the plot parameters for a given medium in scene.plot_eps()."""
global mpl
if mpl is _MPL_NOT_LOADED:
try:
import matplotlib as _mpl
mpl = _mpl
except ImportError:
mpl = None # Triggers grayscale fallback
plot_params = plot_params_structure.copy(update={"linewidth": 0})
if isinstance(medium, AbstractMedium):
if medium.viz_spec is not None:
plot_params = plot_params.override_with_viz_spec(medium.viz_spec)
if alpha is not None:
plot_params = plot_params.copy(update={"alpha": alpha})
if medium.is_pec:
# perfect electrical conductor
plot_params = plot_params.copy(
update={"facecolor": "gold", "edgecolor": "k", "linewidth": 1}
)
elif medium.is_pmc:
# perfect magnetic conductor
plot_params = plot_params.copy(
update={"facecolor": "purple", "edgecolor": "k", "linewidth": 1}
)
elif isinstance(medium, Medium2D):
# 2d material
plot_params = plot_params.copy(update={"edgecolor": "k", "linewidth": 1})
else:
eps_medium = medium._eps_plot(frequency=freq, eps_component=eps_component)
if norm is not None:
color_value = float(norm(eps_medium))
elif eps_min == eps_max:
# Avoid division by zero; use midpoint
color_value = 0.5
else:
color_value = (eps_medium - eps_min) / (eps_max - eps_min)
color_value = min(1.0, max(0.0, color_value))
# Use matplotlib colormap if available, otherwise fallback to grayscale
if mpl is not None:
cmap_name = _get_colormap(reverse=reverse)
cmap = mpl.colormaps[cmap_name]
rgba = tuple(float(component) for component in cmap(color_value))
else:
# Grayscale fallback when matplotlib is unavailable
# Forward (Greys): 0 -> white, 1 -> black
# Reverse (Greys_r): 0 -> black, 1 -> white
grayscale = color_value if reverse else (1.0 - color_value)
rgba = (grayscale, grayscale, grayscale, 1.0)
plot_params = plot_params.copy(update={"facecolor": rgba})
return plot_params
def _plot_shape_structure_eps(
self,
freq: float,
medium: Medium,
shape: Shapely,
eps_min: float,
eps_max: float,
ax: Ax,
reverse: bool = False,
alpha: float | None = None,
eps_component: PermittivityComponent | None = None,
norm: mpl.colors.Normalize | None = None,
) -> Ax:
"""Plot a structure's cross section shape for a given medium, grayscale for permittivity."""
plot_params = self._get_structure_eps_plot_params(
medium=medium,
freq=freq,
eps_min=eps_min,
eps_max=eps_max,
alpha=alpha,
reverse=reverse,
eps_component=eps_component,
norm=norm,
)
ax = self.box.plot_shape(shape=shape, plot_params=plot_params, ax=ax)
return ax
""" Plotting Heat """
[docs]
@equal_aspect
@add_ax_if_none
def plot_heat_charge_property(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
alpha: float | None = None,
cbar: bool = True,
property: str = "heat_conductivity",
ax: Ax = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
) -> Ax:
"""Plot each of scebe's components on a plane defined by one nonzero x,y,z coordinate.
The thermal conductivity is plotted in grayscale based on its value.
Parameters
----------
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
alpha : float = None
Opacity of the structures being plotted.
Defaults to the structure default alpha.
cbar : bool = True
Whether to plot a colorbar for the thermal conductivity.
property : str = "heat_conductivity"
The heat-charge siimulation property to plot. The options are
["heat_conductivity", "electric_conductivity"]
ax : matplotlib.axes._subplots.Axes = None
Matplotlib axes to plot on, if not specified, one is created.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
hlim, vlim = Scene._get_plot_lims(bounds=self.bounds, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
ax = self.plot_structures_heat_charge_property(
cbar=cbar, alpha=alpha, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim, property=property
)
ax = self._set_plot_bounds(bounds=self.bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
return ax
[docs]
@equal_aspect
@add_ax_if_none
def plot_structures_heat_charge_property(
self,
x: float | None = None,
y: float | None = None,
z: float | None = None,
alpha: float | None = None,
cbar: bool = True,
property: str = "heat_conductivity",
reverse: bool = False,
ax: Ax = None,
hlim: tuple[float, float] | None = None,
vlim: tuple[float, float] | None = None,
) -> Ax:
"""Plot each of scene's structures on a plane defined by one nonzero x,y,z coordinate.
The thermal conductivity is plotted in grayscale based on its value.
Parameters
----------
x : float = None
position of plane in x direction, only one of x, y, z must be specified to define plane.
y : float = None
position of plane in y direction, only one of x, y, z must be specified to define plane.
z : float = None
position of plane in z direction, only one of x, y, z must be specified to define plane.
reverse : bool = False
If ``False``, the highest permittivity is plotted in black.
If ``True``, it is plotteed in white (suitable for black backgrounds).
cbar : bool = True
Whether to plot a colorbar for the relative permittivity.
alpha : float = None
Opacity of the structures being plotted.
Defaults to the structure default alpha.
ax : matplotlib.axes._subplots.Axes = None
Matplotlib axes to plot on, if not specified, one is created.
hlim : tuple[float, float] = None
The x range if plotting on xy or xz planes, y range if plotting on yz plane.
vlim : tuple[float, float] = None
The z range if plotting on xz or yz planes, y plane if plotting on xy plane.
Returns
-------
matplotlib.axes._subplots.Axes
The supplied or created matplotlib axes.
"""
structures = self.sorted_structures
# alpha is None just means plot without any transparency
if alpha is None:
alpha = 1
if alpha <= 0:
return ax
if alpha < 1:
axis, position = Box.parse_xyz_kwargs(x=x, y=y, z=z)
center = Box.unpop_axis(position, (0, 0), axis=axis)
size = Box.unpop_axis(0, (inf, inf), axis=axis)
plane = Box(center=center, size=size)
medium_shapes = self._filter_structures_plane_medium(
structures=structures, plane=plane, section_tolerance_2d=True
)
else:
structures = [self.background_structure, *list(structures)]
medium_shapes = self._get_structures_2dbox(
structures=structures, x=x, y=y, z=z, hlim=hlim, vlim=vlim
)
property_val_min, property_val_max = self.heat_charge_property_bounds(property=property)
for medium, shape in medium_shapes:
ax = self._plot_shape_structure_heat_charge_property(
alpha=alpha,
medium=medium,
property_val_min=property_val_min,
property_val_max=property_val_max,
reverse=reverse,
shape=shape,
ax=ax,
property=property,
)
if cbar:
label = ""
if property == "heat_conductivity":
label = f"Thermal conductivity ({THERMAL_CONDUCTIVITY})"
elif property == "electric_conductivity":
label = f"Electric conductivity ({CONDUCTIVITY})"
self._add_cbar(
vmin=property_val_min,
vmax=property_val_max,
label=label,
cmap=STRUCTURE_HEAT_COND_CMAP,
ax=ax,
)
# clean up the axis display
axis, _ = Box.parse_xyz_kwargs(x=x, y=y, z=z)
ax = self.box.add_ax_lims(axis=axis, ax=ax)
ax = self._set_plot_bounds(bounds=self.bounds, ax=ax, x=x, y=y, z=z, hlim=hlim, vlim=vlim)
# Add the default axis labels, tick labels, and title
ax = Box.add_ax_labels_and_title(
ax=ax, x=x, y=y, z=z, plot_length_units=self.plot_length_units
)
return ax
[docs]
def heat_charge_property_bounds(self, property: str) -> tuple[float, float]:
"""Compute range of the heat-charge simulation property present in the scene.
Returns
-------
tuple[float, float]
Minimal and maximal values of thermal conductivity in scene.
"""
medium_list = [self.medium, *list(self.mediums)]
if property == "heat_conductivity":
SolidType = (SolidSpec, SolidMedium)
medium_list = [
medium for medium in medium_list if isinstance(medium.heat_spec, SolidType)
]
cond_list = [
_heat_conductivity_scalar(medium.heat_spec.conductivity) for medium in medium_list
]
elif property == "electric_conductivity":
cond_mediums = [
medium for medium in medium_list if isinstance(medium.charge, ChargeConductorMedium)
]
cond_list = [medium.charge.conductivity for medium in cond_mediums]
if len(cond_list) == 0:
cond_list = [0]
cond_min = min(cond_list)
cond_max = max(cond_list)
return cond_min, cond_max
def _get_structure_heat_charge_property_plot_params(
self,
medium: Medium,
property_val_min: float,
property_val_max: float,
reverse: bool = False,
alpha: float | None = None,
property: str = "heat_conductivity",
) -> PlotParams:
"""Constructs the plot parameters for a given medium in
scene.plot_heat_charge_property().
"""
plot_params = plot_params_structure.copy(update={"linewidth": 0})
if hasattr(medium, "viz_spec"):
if medium.viz_spec is not None:
plot_params = plot_params.override_with_viz_spec(medium.viz_spec)
if alpha is not None:
plot_params = plot_params.copy(update={"alpha": alpha})
cond_medium = None
SolidType = (SolidSpec, SolidMedium)
if property == "heat_conductivity" and isinstance(medium.heat_spec, SolidType):
cond_medium = _heat_conductivity_scalar(medium.heat_spec.conductivity)
elif property == "electric_conductivity" and isinstance(
medium.charge, ChargeConductorMedium
):
cond_medium = medium.charge.conductivity
elif property == "doping":
cond_medium = None
if cond_medium is not None:
delta_cond = cond_medium - property_val_min
delta_cond_max = property_val_max - property_val_min + 1e-5 * property_val_min
cond_fraction = delta_cond / delta_cond_max
color = cond_fraction if reverse else 1 - cond_fraction
plot_params = plot_params.copy(update={"facecolor": str(color)})
else:
plot_params = plot_params_fluid
if alpha is not None:
plot_params = plot_params.copy(update={"alpha": alpha})
return plot_params
def _plot_shape_structure_heat_charge_property(
self,
medium: Medium,
shape: Shapely,
property_val_min: float,
property_val_max: float,
property: str,
ax: Ax,
reverse: bool = False,
alpha: float | None = None,
) -> Ax:
"""Plot a structure's cross section shape for a given medium, grayscale for thermal
conductivity.
"""
plot_params = self._get_structure_heat_charge_property_plot_params(
medium=medium,
property_val_min=property_val_min,
property_val_max=property_val_max,
alpha=alpha,
reverse=reverse,
property=property,
)
ax = self.box.plot_shape(shape=shape, plot_params=plot_params, ax=ax)
return ax
""" Misc """
[docs]
def perturbed_mediums_copy(
self,
temperature: CustomSpatialDataType | None = None,
electron_density: CustomSpatialDataType | None = None,
hole_density: CustomSpatialDataType | None = None,
interp_method: InterpMethod = "linear",
) -> Self:
"""Return a copy of the scene with heat and/or charge data applied to all mediums
that have perturbation models specified. That is, such mediums will be replaced with
spatially dependent custom mediums that reflect perturbation effects. Any of temperature,
electron_density, and hole_density can be ``None``. All provided fields must have identical
coords.
Note
----
The resulting custom mediums are given unique names based on the structure they belong to.
If the original medium has a name (e.g., ``"Si"``), the new medium name will be
``"Si[structure_name]"`` if the structure has a name, or ``"Si[structures[index]]"``
otherwise. This ensures unique medium names across the scene.
Parameters
----------
temperature : Union[:class:`.SpatialDataArray`, :class:`.TriangularGridDataset`, :class:`.TetrahedralGridDataset`] = None
Temperature field data.
electron_density : Union[:class:`.SpatialDataArray`, :class:`.TriangularGridDataset`, :class:`.TetrahedralGridDataset`] = None
Electron density field data.
hole_density : Union[:class:`.SpatialDataArray`, :class:`.TriangularGridDataset`, :class:`.TetrahedralGridDataset`] = None
Hole density field data.
interp_method : :class:`.InterpMethod`, optional
Interpolation method to obtain heat and/or charge values that are not supplied
at the Yee grids.
Returns
-------
:class:`.Scene`
Simulation after application of heat and/or charge data.
"""
scene_dict = self.model_dump()
structures = self.sorted_structures
array_dict = {
"temperature": temperature,
"electron_density": electron_density,
"hole_density": hole_density,
}
# For each structure made of mediums with perturbation models, convert those mediums into
# spatially dependent mediums by selecting minimal amount of heat and charge data points
# covering the structure, and create a new structure containing the resulting custom medium
new_structures = []
for s_ind, structure in enumerate(structures):
med = structure.medium
if isinstance(med, AbstractPerturbationMedium):
# get structure's bounding box
bounds = structure.geometry.bounds
# for each structure select a minimal subset of data that covers it
restricted_arrays = {}
for name, array in array_dict.items():
if array is not None:
restricted_arrays[name] = array.sel_inside(bounds)
# check provided data fully cover structure
if not array.does_cover(bounds):
log.warning(
f"Provided '{name}' does not fully cover structures[{s_ind}]."
)
new_medium = med.perturbed_copy(**restricted_arrays, interp_method=interp_method)
# Generate unique medium name based on structure to avoid duplicate name warnings.
# Only rename if a new medium was actually created (not when perturbed_copy returns
# the original unchanged due to no perturbation data being provided).
if new_medium is not med and new_medium.name is not None:
suffix = structure.name if structure.name else f"structures[{s_ind}]"
new_medium = new_medium.updated_copy(name=f"{new_medium.name}[{suffix}]")
new_structure = structure.updated_copy(medium=new_medium)
new_structures.append(new_structure)
else:
new_structures.append(structure)
scene_dict["structures"] = new_structures
# do the same for background medium if it a medium with perturbation models.
med = self.medium
if isinstance(med, AbstractPerturbationMedium):
scene_dict["medium"] = med.perturbed_copy(**array_dict, interp_method=interp_method)
return Scene.model_validate(scene_dict)
[docs]
def doping_bounds(self) -> tuple[list[float], list[float]]:
"""Get the maximum and minimum of the doping"""
acceptors_lims = [np.inf, -np.inf]
donors_lims = [np.inf, -np.inf]
for struct in self.all_structures:
if isinstance(struct.medium.charge, SemiconductorMedium):
electric_spec = struct.medium.charge
for doping, limits in zip(
[electric_spec.N_a, electric_spec.N_d], [acceptors_lims, donors_lims]
):
if isinstance(doping, float):
if doping < limits[0]:
limits[0] = doping
if doping > limits[1]:
limits[1] = doping
# NOTE: This will be deprecated.
if isinstance(doping, SpatialDataArray):
min_value = np.min(doping.data.flatten())
max_value = np.max(doping.data.flatten())
if min_value < limits[0]:
limits[0] = min_value
if max_value > limits[1]:
limits[1] = max_value
if isinstance(doping, tuple):
for doping_box in doping:
if isinstance(doping_box, ConstantDoping):
if doping_box.concentration < limits[0]:
limits[0] = doping_box.concentration
if doping_box.concentration > limits[1]:
limits[1] = doping_box.concentration
if isinstance(doping_box, GaussianDoping):
if doping_box.ref_con < limits[0]:
limits[0] = doping_box.ref_con
if doping_box.concentration > limits[1]:
limits[1] = doping_box.concentration
if isinstance(doping_box, CustomDoping):
min_value = np.min(doping_box.concentration.data.flatten())
max_value = np.max(doping_box.concentration.data.flatten())
if min_value < limits[0]:
limits[0] = min_value
if max_value > limits[1]:
limits[1] = max_value
# make sure we have recorded some values. Otherwise, set to 0
if np.isinf(acceptors_lims[0]):
acceptors_lims[0] = 0
if np.isinf(acceptors_lims[1]):
acceptors_lims[1] = 0
if np.isinf(donors_lims[0]):
donors_lims[0] = 0
if np.isinf(donors_lims[1]):
donors_lims[1] = 0
return acceptors_lims, donors_lims
[docs]
def doping_absolute_minimum(self) -> tuple[float, float]:
"""Get the absolute minimum values of the doping concentrations.
Returns
-------
Tuple[float, float]
Absolute minimum values for acceptors and donors respectively.
"""
# Use more reasonable initial values
acceptors_abs_min = np.inf
donors_abs_min = np.inf
for struct in self.all_structures:
if isinstance(struct.medium.charge, SemiconductorMedium):
electric_spec = struct.medium.charge
# Process acceptors
acceptors_min = self._get_absolute_minimum_from_doping(electric_spec.N_a)
if acceptors_min < acceptors_abs_min:
acceptors_abs_min = acceptors_min
# Process donors
donors_min = self._get_absolute_minimum_from_doping(electric_spec.N_d)
if donors_min < donors_abs_min:
donors_abs_min = donors_min
return acceptors_abs_min, donors_abs_min
def _get_absolute_minimum_from_doping(
self, doping: float | SpatialDataArray | tuple[DopingBoxType, ...]
) -> float:
"""Helper method to get absolute minimum from a single doping specification.
Parameters
----------
doping : Union[float, SpatialDataArray, tuple]
Doping specification to analyze.
Returns
-------
float
Absolute minimum value found in the doping specification.
"""
if isinstance(doping, float):
return np.abs(doping)
# NOTE: This will be deprecated.
if isinstance(doping, SpatialDataArray):
return np.min(np.abs(doping.data.flatten()))
if isinstance(doping, tuple):
min_values = []
for doping_box in doping:
if isinstance(doping_box, ConstantDoping):
min_values.append(np.abs(doping_box.concentration))
elif isinstance(doping_box, GaussianDoping):
min_values.append(np.abs(doping_box.ref_con))
elif isinstance(doping_box, CustomDoping):
min_values.append(np.min(np.abs(doping_box.concentration.data.flatten())))
return min(min_values) if min_values else np.inf
return np.inf
def _pcolormesh_shape_doping_box(
self,
x: float,
y: float,
z: float,
alpha: float,
medium: Medium,
shape: Shapely,
ax: Ax,
plt_type: str = "doping",
norm: mpl.colors.Normalize | None = None,
) -> None:
"""
Plot shape made of structure defined with doping.
plt_type accepts ["doping", "N_a", "N_d"]
"""
coords = "xyz"
normal_axis_ind, normal_position = Box.parse_xyz_kwargs(x=x, y=y, z=z)
_normal_axis, _plane_axes = Box.pop_axis(coords, normal_axis_ind)
# make grid for eps interpolation
# we will do this by combining shape bounds and points where custom eps is provided
shape_bounds = shape.bounds
rmin, rmax = [*shape_bounds[:2]], [*shape_bounds[2:]]
rmin.insert(normal_axis_ind, normal_position)
rmax.insert(normal_axis_ind, normal_position)
# for the time being let's assume we'll always need to generate a mesh
plane_axes_inds = [0, 1, 2]
plane_axes_inds.pop(normal_axis_ind)
# build grid
N = 100
coords_2D = [np.linspace(rmin[d], rmax[d], N) for d in plane_axes_inds]
X, Y = np.meshgrid(coords_2D[0], coords_2D[1], indexing="ij")
struct_doping = [
np.zeros(X.shape), # let's use 0 for N_a
np.zeros(X.shape), # and 1 for N_d
]
electric_spec = medium.charge
for n, doping in enumerate([electric_spec.N_a, electric_spec.N_d]):
if isinstance(doping, float):
struct_doping[n] = struct_doping[n] + doping
# NOTE: This will be deprecated.
if isinstance(doping, SpatialDataArray):
struct_coords = {"xyz"[d]: coords_2D[i] for i, d in enumerate(plane_axes_inds)}
data_2D = doping
# check whether the provided doping data is 2 or 3D
data_is_2d = any(dim_size <= 1 for _, dim_size in doping.sizes.items())
if not data_is_2d:
selector = {"xyz"[normal_axis_ind]: normal_position}
data_2D = doping.sel(**selector)
contrib = data_2D.interp(
**struct_coords,
method="nearest",
kwargs={"bounds_error": False, "fill_value": 0},
)
struct_doping[n] = struct_doping[n] + contrib
# Handle doping boxes
if isinstance(doping, tuple):
boxes = [b for b in doping if isinstance(b, DopingBoxType.__args__)]
if boxes:
# Use the same seam-aware aggregation as the solver so the
# rendered field matches the solved one (interiors SUM, abutting
# seams dedupe via last-in-list-wins, 1-ULP misses recovered).
flat = [None, None, None]
flat[plane_axes_inds[0]] = X.ravel()
flat[plane_axes_inds[1]] = Y.ravel()
flat[normal_axis_ind] = np.full(X.size, normal_position)
contrib = aggregate_doping_seam_aware(boxes, *flat)
struct_doping[n] = struct_doping[n] + contrib.reshape(X.shape)
if plt_type == "doping":
struct_doping_to_plot = struct_doping[0] - struct_doping[1]
elif plt_type == "N_a":
struct_doping_to_plot = struct_doping[0]
elif plt_type == "N_d":
struct_doping_to_plot = struct_doping[1]
ax.pcolormesh(
X,
Y,
struct_doping_to_plot,
clip_path=(polygon_path(shape), ax.transData),
cmap=HEAT_SOURCE_CMAP,
alpha=alpha,
clip_box=ax.bbox,
norm=norm,
)
[docs]
def plot_3d(self, width: int = 800, height: int = 800) -> None:
"""Render 3D plot of ``Scene`` (in jupyter notebook only).
Parameters
----------
width : float = 800
width of the 3d view dom's size
height : float = 800
height of the 3d view dom's size
"""
return plot_scene_3d(self, width=width, height=height)
