feat: Added the Bloch boundary condition.

A very healthy amount of research on how to choose the Bloch vector was
performed.
It is encapsulated not only in the documentation, but also in the modes
available for how to derive one that fits a given simulation.

The theory of the Bloch boundaries can really bite you, and the hope is
that by focusing on such invalid-usage-prevention, a lot of time can be
saved in the sim design stages.

TODO.md

@@ -1,12 +1,5 @@
 # Working TODO
 - [x] Wave Constant
-- Bounds
-	- [x] Boundary Conds
-	- [x] PML
-	- [x] PEC
-	- [x] PMC
-	- [ ] Bloch
-	- [ ] Absorbing
 - Sources
 	- [ ] Temporal Shapes / Continuous Wave Temporal Shape
 	- [ ] Temporal Shapes / Symbolic Temporal Shape
@@ -200,9 +193,11 @@
 - [x] Boundary Conds
 - [x] Boundary Cond / PML Bound Cond
 	- [ ] 1D plot visualizing the effect of parameters on a 1D wave function
-- [ ] Boundary Cond / Bloch Bound Cond
-	- [ ] Implement "simple" mode aka "periodic" mode in Tidy3D
-- [ ] Boundary Cond / Absorbing Bound Cond
+- [x] Boundary Cond / Bloch Bound Cond
+	- [x] Implement "simple" mode aka "periodic" mode in Tidy3D
+	- [ ] 1D plot visualizing the effect of parameters on a 1D wave function
+- [x] Boundary Cond / Absorbing Bound Cond
+	- [ ] 1D plot visualizing the effect of parameters on a 1D wave function
 
 ## Monitors
 - [x] EH Field Monitor

src/blender_maxwell/node_trees/maxwell_sim_nodes/contracts/__init__.py

@@ -1,25 +1,25 @@
 from blender_maxwell.contracts import (
-	BLClass,
-	BLColorRGBA,
-	BLEnumElement,
-	BLEnumID,
-	BLIcon,
-	BLIconSet,
-	BLIDStruct,
-	BLKeymapItem,
-	BLModifierType,
-	BLNodeTreeInterfaceID,
-	BLOperatorStatus,
-	BLPropFlag,
-	BLRegionType,
-	BLSpaceType,
-	KeymapItemDef,
-	ManagedObjName,
-	OperatorType,
-	PanelType,
-	PresetName,
-	SocketName,
-	addon,
+    BLClass,
+    BLColorRGBA,
+    BLEnumElement,
+    BLEnumID,
+    BLIcon,
+    BLIconSet,
+    BLIDStruct,
+    BLKeymapItem,
+    BLModifierType,
+    BLNodeTreeInterfaceID,
+    BLOperatorStatus,
+    BLPropFlag,
+    BLRegionType,
+    BLSpaceType,
+    KeymapItemDef,
+    ManagedObjName,
+    OperatorType,
+    PanelType,
+    PresetName,
+    SocketName,
+    addon,
 )
 
 from .bl_socket_types import BLSocketInfo, BLSocketType
@@ -27,20 +27,20 @@ from .category_labels import NODE_CAT_LABELS
 from .category_types import NodeCategory
 from .flow_events import FlowEvent
 from .flow_kinds import (
-	ArrayFlow,
-	CapabilitiesFlow,
-	FlowKind,
-	InfoFlow,
-	LazyArrayRangeFlow,
-	LazyValueFuncFlow,
-	ParamsFlow,
-	ValueFlow,
+    ArrayFlow,
+    CapabilitiesFlow,
+    FlowKind,
+    InfoFlow,
+    LazyArrayRangeFlow,
+    LazyValueFuncFlow,
+    ParamsFlow,
+    ValueFlow,
 )
 from .flow_signals import FlowSignal
 from .icons import Icon
 from .mobj_types import ManagedObjType
 from .node_types import NodeType
-from .sim_types import BoundCondType
+from .sim_types import BoundCondType, SimSpaceAxis
 from .socket_colors import SOCKET_COLORS
 from .socket_types import SocketType
 from .tree_types import TreeType
@@ -79,6 +79,7 @@ __all__ = [
 	'BLSocketType',
 	'NodeType',
 	'BoundCondType',
+	'SimSpaceAxis',
 	'NodeCategory',
 	'NODE_CAT_LABELS',
 	'ManagedObjType',

src/blender_maxwell/node_trees/maxwell_sim_nodes/contracts/flow_kinds.py

@@ -13,9 +13,12 @@ import sympy as sp
 import sympy.physics.units as spu
 
 from blender_maxwell.utils import extra_sympy_units as spux
+from blender_maxwell.utils import logger
 
 from .socket_types import SocketType
 
+log = logger.get(__name__)
+
 
 class FlowKind(enum.StrEnum):
 	"""Defines a kind of data that can flow between nodes.
@@ -87,17 +90,28 @@ class CapabilitiesFlow:
 	active_kind: FlowKind
 
 	is_universal: bool = False
+
+	# == Constraint
 	must_match: dict[str, typ.Any] = dataclasses.field(default_factory=dict)
 
+	# ∀b (b ∈ A) Constraint
+	## A: allow_any
+	## b∈B: present_any
+	allow_any: set[typ.Any] = dataclasses.field(default_factory=set)
+	present_any: set[typ.Any] = dataclasses.field(default_factory=set)
+
 	def is_compatible_with(self, other: typ.Self) -> bool:
 		return other.is_universal or (
 			self.socket_type == other.socket_type
 			and self.active_kind == other.active_kind
+			# == Constraint
 			and all(
 				name in other.must_match
 				and self.must_match[name] == other.must_match[name]
 				for name in self.must_match
 			)
+			# ∀b (b ∈ A) Constraint
+			and self.present_any.issubset(other.allow_any)
 		)
 
 

src/blender_maxwell/node_trees/maxwell_sim_nodes/contracts/sim_types.py

@@ -6,6 +6,55 @@ import typing as typ
 import tidy3d as td
 
 
+## TODO: Sim Domain type, w/pydantic checks!
+
+class SimSpaceAxis(enum.StrEnum):
+	"""The axis labels of the global simulation coordinate system."""
+
+	X = enum.auto()
+	Y = enum.auto()
+	Z = enum.auto()
+
+	@staticmethod
+	def to_name(v: typ.Self) -> str:
+		"""Convert the enum value to a human-friendly name.
+
+		Notes:
+			Used to print names in `EnumProperty`s based on this enum.
+
+		Returns:
+			A human-friendly name corresponding to the enum value.
+		"""
+		SSA = SimSpaceAxis
+		return {
+			SSA.X: 'x',
+			SSA.Y: 'y',
+			SSA.Z: 'z',
+		}[v]
+
+	@staticmethod
+	def to_icon(_: typ.Self) -> str:
+		"""Convert the enum value to a Blender icon.
+
+		Notes:
+			Used to print icons in `EnumProperty`s based on this enum.
+
+		Returns:
+			A human-friendly name corresponding to the enum value.
+		"""
+		return ''
+
+	@property
+	def axis(self) -> int:
+		"""Deduce the integer index of the axis.
+
+		Returns:
+			The integer index of the axis.
+		"""
+		SSA = SimSpaceAxis
+		return {SSA.X: 0, SSA.Y: 1, SSA.Z: 2}[self]
+
+
 class BoundCondType(enum.StrEnum):
 	r"""A type of boundary condition, applied to a half-axis of a simulation domain.
 

src/blender_maxwell/node_trees/maxwell_sim_nodes/nodes/bounds/bound_cond_nodes/__init__.py

@@ -1,16 +1,16 @@
 from . import (
 	absorbing_bound_cond,
-	# bloch_bound_cond,
+	bloch_bound_cond,
 	pml_bound_cond,
 )
 
 BL_REGISTER = [
 	*pml_bound_cond.BL_REGISTER,
-	# *bloch_bound_cond.BL_REGISTER,
+	*bloch_bound_cond.BL_REGISTER,
 	*absorbing_bound_cond.BL_REGISTER,
 ]
 BL_NODES = {
 	**pml_bound_cond.BL_NODES,
-	# **bloch_bound_cond.BL_NODES,
+	**bloch_bound_cond.BL_NODES,
 	**absorbing_bound_cond.BL_NODES,
 }

src/blender_maxwell/node_trees/maxwell_sim_nodes/nodes/bounds/bound_cond_nodes/absorbing_bound_cond.py

@@ -145,5 +145,5 @@ BL_REGISTER = [
 	AdiabAbsorbBoundCondNode,
 ]
 BL_NODES = {
-	ct.NodeType.AdiabAbsorbBoundCond: (ct.NodeCategory.MAXWELLSIM_BOUNDS_BOUNDCONDS)
+	ct.NodeType.AdiabAbsorbBoundCond: (ct.NodeCategory.MAXWELLSIM_BOUNDS)
 }

src/blender_maxwell/node_trees/maxwell_sim_nodes/nodes/bounds/bound_cond_nodes/bloch_bound_cond.py

@@ -1,5 +1,236 @@
+"""Implements `BlochBoundCondNode`."""
+
+import typing as typ
+
+import bpy
+import tidy3d as td
+
+from blender_maxwell.utils import bl_cache, logger
+
+from .... import contracts as ct
+from .... import sockets
+from ... import base, events
+
+log = logger.get(__name__)
+
+
+class BlochBoundCondNode(base.MaxwellSimNode):
+	r"""A boundary condition that declares an "infinitely repeating" window, by applying Bloch's theorem to accurately describe how a boundary would behave if it were interacting with an infinitely repeating simulation structure.
+
+	# Theory
+	In the simplest case, aka. a normal-incident plane wave, the symmetries of electromagnetic wave propagation behave exactly as expected: Copy-paste the wavevector, but at opposite corners, as part of the FDTD neighbor-cell-update.
+	The moment this plane wave becomes angled, however, this "naive" method will cause **the phase of the periodically propagated fields to diverge from reality**.
+
+	With a bit of hand-waving, this is a natural thing: Fundamentally, the distance from each point on an angled plane wave to the boundary must vary, and if the phase is distance-dependent, then the phase must vary across the boundary.
+
+	Unfortunately, all of the explicitly-defined ways of describing how exactly to correct for this phenomenon depend on not only on what is being simulated, but on what is being studied.
+	The good news is, there are options.
+
+	## A Bloch of a Thing
+	A physicist named Felix Bloch came up with a theorem to help constrain how "wave-like things in periodic stuff" can be thought about, and it looks like
+
+	$$
+		\psi(\mathbf{r}) = u(\mathbf{r}) \cdot \exp_{\mathbb{C}}(\mathbf{k} \cdot \mathbf{r})
+	$$
+
+	for:
+
+	- $\psi$: A wave function (in general, satisfying the Schrödinger equation, but in this context, satisfying Maxwell's equations)
+	- $\mathbf{r}$: A position in 3D space.
+	- $\u$: Some periodic function mapping 3D space to a value. In this context, this might be a 3D function representing our simulation structures.
+	- $\mathbf{k}$: The "Bloch vector", of which there is guaranteed to be at least one, **but of which there may be many**.
+
+	At this point, it becomes interesting to note that pretty much _everything_ is, in fact, a "wave-like thing", so long as "the thing" is small enough.
+	Many such "periodically structured things", which form entire fields of study, can indeed be modelled using this single function:
+
+	- **Photonic Crystals**: The optical properties of many materials can be quite concisely encapsulated by placing regularly placed structures (of sub-wavelength size) within lattice-like structures.
+	- **Phononic Crystals**: A class of metamaterial that can be parameterized and optimized for its acoustic properties, purely by analyzing its periodic behavior, with applications ranging from interesting acoustic devices to seismic modelling.
+
+	## Modes of an Excited Nature
+	For a choice of $u$ (representing the simulation structure), there may be _many continuous_ choices of $\mathbf{k}$ that satisfy the Bloch theorem.
+	Similarly, for a particular choice of $\mathbf{k}$, there may be _several discrete_ particular solutions of the given wave function.
+
+	Thus, we come full circle: **Fully encapsulating** the wave-interactions of a periodic structure requires knowing its behavior at **all valid wave vectors**.
+	It is a sort of deeper truth, that any particular simulation of a unit cell cannot elicit the full story of _how a structure behaves_, since a particular choice of $\mathbf{k}$ must always inevitably be made as part of defining the simulation.
+
+	## Designing Periodically
+	With this insight in mind, we can now design simulations of periodic structures that properly account for the modalities imposed by particular $\mathbf{k}$ choices:
+
+	- **Only Rely on Real Fields**: If only the real parts of the fields are important, then the choice of $\mathbf{k}$ might not matter.
+		Remember, the symptom of needing to understand $\mathbf{k}$ is the phase-shift; if the phase-shift does not matter, then altering the Bloch vector won't change anything.
+		**Be careful**, though, and make sure to validate that the Bloch vector truly doesn't change anything.
+	- **Normal-Injected Plane Waves**: If fields generally only propagate in the normal direction, then again, choices of $\mathbf{k}$ might not matter.
+		Again, phase-shifting due to periodic behavior mainly happens when propagation occurs at grazing angles.
+		Again, **be careful**, and make sure to validate that ex. the Poynting vector truly isn't hitting the boundaries at too-grazing angles.
+	- **Angularly Injected Plane Waves**: If the injected plane wave is known, then we can directly compute a reasonable Bloch vector from the angle and boundary-axis-projected size of the plane wave source.
+		This selection of $\mathbf{k}$
+	- **Brute-Force Bloch-Vector Sweep**: If the nature of a periodic structure needs to be uncovered, and there's no special constraints to rely on, then it would be rightfully tempting to just sweep over all $\mathbf{k}$s, and run a complete simluation for each.
+		By going a step further, and plotting the energy of resonance frequencies noticed at each wave vector (just place point dipoles at random), one might stumble into a "band diagram" describing the possible energy states of electrons at each wave vector.
+
+	In general, these form a very sensible starting point for how to select Bloch vectors for productive use in the simulation.
+
+	NOTE: The Bloch vector is generally represented not as a vector, but as a single phase-shift per boundary axis unit length, mainly for convenience.
+
+	## Further Reading
+	- <https://optics.ansys.com/hc/en-us/articles/360041566614-Rectangular-Photonic-Crystal-Bandstructure>
+	- <https://docs.flexcompute.com/projects/tidy3d/en/v2.1.0/notebooks/Bandstructure.html>
+	- <https://en.wikipedia.org/wiki/Electronic_band_structure>
+	- <https://en.wikipedia.org/wiki/Brillouin_zone>
+	- <https://en.wikipedia.org/wiki/Bloch%27s_theorem>
+
+	Notes:
+		In the naive case, it is presumed that the choice of Bloch vector doesn't matter; therefore it is set to 0.
+
+	Socket Sets:
+		Naive: Specify a Bloch boundary condition where phase shift doesn't matter, and is thus set such that no phase-shift occurs.
+			This is the simplest (and cheapest) mechanism, which merely copy-pastes propagating waves at opposing sides of the simulation.
+			However, **this should not be used for angled plane waves**, as the phase-shift of a propagating angled plane wave **will be wrong**.
+		Source-Derived: Derive a Bloch vector that will be generally correct for a directed source, within a particular choice of axis on a particular simulation domain.
+			**Phase shift correctness is only guaranteed valid for the central frequency of the source**.
+			Thus, a narrow-band source is strongly recommended.
+		Bloch Vector: Specify a true Bloch boundary condition, including the **phase shift per unit length** (aka. the magnitude of the Bloch vector).
+			While the most flexible, **the appropriate choice for this value source of this value depends entirely on what is being simulated**.
+	"""
+
+	node_type = ct.NodeType.BlochBoundCond
+	bl_label = 'Bloch Bound Cond'
+
+	####################
+	# - Sockets
+	####################
+	input_socket_sets: typ.ClassVar = {
+		'Naive': {},
+		'Source-Derived': {
+			'Angled Source': sockets.MaxwellSourceSocketDef(),
+			## TODO: Constrain to gaussian beam, plane wafe, and tfsf
+			'Sim Domain': sockets.MaxwellSimDomainSocketDef(),
+		},
+		'Manual': {
+			'Bloch Vector': sockets.ExprSocketDef(),
+		},
+	}
+	output_sockets: typ.ClassVar = {
+		'BC': sockets.MaxwellBoundCondSocketDef(),
+	}
+
+	####################
+	# - Properties
+	####################
+	valid_sim_axis: ct.SimSpaceAxis = bl_cache.BLField(ct.SimSpaceAxis.X, prop_ui=True)
+
+	####################
+	# - UI
+	####################
+	def draw_props(self, _: bpy.types.Context, layout: bpy.types.UILayout) -> None:
+		if self.active_socket_set == 'Source-Derived':
+			layout.prop(self, self.blfields['valid_sim_axis'], expand=True)
+
+	def draw_info(self, _: bpy.types.Context, layout: bpy.types.UILayout) -> None:
+		if self.active_socket_set == 'Manual':
+			box = layout.box()
+			row = box.row()
+			row.alignment = 'CENTER'
+			row.label(text='Interpretation')
+
+			# Split
+			split = box.split(factor=0.6, align=False)
+
+			## LHS: Parameter Names
+			col = split.column()
+			col.alignment = 'RIGHT'
+			col.label(text='Bloch Vec:')
+
+			## RHS: Parameter Units
+			col = split.column()
+			col.label(text='2π/Δℓ')
+
+	####################
+	# - Events
+	####################
+	@events.on_value_changed(
+		prop_name={'active_socket_set', 'valid_sim_axis'},
+		run_on_init=True,
+		props={'active_socket_set', 'valid_sim_axis'},
+	)
+	def on_valid_sim_axis_changed(self, props):
+		"""For the source-derived socket set, synchronized the output socket's axis compatibility with the axis onto which the Bloch vector is computed.
+
+		The net result should be that invalid use of the Bloch boundary condition in a particular axis should be rejected.
+
+		- **Source-Derived**: Since the Bloch vector is computed between the source and the axis that this boundary is applied to, the output socket must be altered to **only** declare compatibility with that axis.
+		- **`*`**: Normalize the output socket axis validity to ensure that the boundary condition can be applied to any axis.
+		"""
+		if props['active_socket_set'] == 'Source-Derived':
+			self.outputs['BC'].present_axes = {props['valid_sim_axis']}
+			self.outputs['BC'].remove_invalidated_links()
+		else:
+			self.outputs['BC'].present_axes = {
+				ct.SimSpaceAxis.X,
+				ct.SimSpaceAxis.Y,
+				ct.SimSpaceAxis.Z,
+			}
+
+	####################
+	# - Output
+	####################
+	@events.computes_output_socket(
+		'BC',
+		props={'active_socket_set', 'valid_sim_axis'},
+		input_sockets={
+			'Angled Source',
+			'Sim Domain',
+			'Bloch Vector',
+		},
+		input_sockets_optional={
+			'Angled Source': True,
+			'Sim Domain': True,
+			'Bloch Vector': True,
+		},
+	)
+	def compute_bloch_bound_cond(
+		self, props, input_sockets
+	) -> td.Periodic | td.BlochBoundary:
+		r"""Computes the Bloch boundary condition.
+
+		- **Naive**: Set the Bloch vector to 0 by returning a `td.Periodic`.
+		- **Source-Derived**: Derive the Bloch vector from the source, simulation domain, and choice of axis.
+			The Bloch boundary axis **must** be orthogonal to the source's injection axis.
+		- **Manual**: Set the Bloch vector to the user-specified value.
+		"""
+		log.debug(
+			'%s: Computing Bloch Boundary Condition (Socket Set = %s)',
+			self.sim_node_name,
+			props['active_socket_set'],
+		)
+
+		# Naive
+		if props['active_socket_set'] == 'Naive':
+			return td.Periodic()
+
+		# Source-Derived
+		if props['active_socket_set'] == 'Naive':
+			sim_domain = input_sockets['Sim Domain']
+			valid_sim_axis = props['valid_sim_axis']
+
+			has_sim_domain = not ct.FlowSignal.check(sim_domain)
+
+			if has_sim_domain:
+				return td.BlochBoundary.from_source(
+					source=input_sockets['Angled Source'],
+					domain_size=sim_domain['size'][valid_sim_axis.axis],
+					axis=valid_sim_axis.axis,
+					medium=sim_domain['medium'],
+				)
+			return ct.FlowSignal.FlowPending
+
+		# Manual
+		return td.BlochBoundary(bloch_vec=input_sockets['Bloch Vector'])
+
+
 ####################
 # - Blender Registration
 ####################
-BL_REGISTER = []
-BL_NODES = {}
+BL_REGISTER = [
+	BlochBoundCondNode,
+]
+BL_NODES = {ct.NodeType.BlochBoundCond: (ct.NodeCategory.MAXWELLSIM_BOUNDS)}

src/blender_maxwell/node_trees/maxwell_sim_nodes/nodes/bounds/bound_cond_nodes/pml_bound_cond.py

@@ -187,4 +187,4 @@ class PMLBoundCondNode(base.MaxwellSimNode):
 BL_REGISTER = [
 	PMLBoundCondNode,
 ]
-BL_NODES = {ct.NodeType.PMLBoundCond: (ct.NodeCategory.MAXWELLSIM_BOUNDS_BOUNDCONDS)}
+BL_NODES = {ct.NodeType.PMLBoundCond: (ct.NodeCategory.MAXWELLSIM_BOUNDS)}

src/blender_maxwell/node_trees/maxwell_sim_nodes/nodes/bounds/bound_conds.py

@@ -13,6 +13,9 @@ from .. import base, events
 log = logger.get(__name__)
 
 
+SSA = ct.SimSpaceAxis
+
+
 class BoundCondsNode(base.MaxwellSimNode):
 	"""Provides a hub for joining custom simulation domain boundary conditions by-axis."""
 
@@ -24,49 +27,49 @@ class BoundCondsNode(base.MaxwellSimNode):
 	####################
 	input_socket_sets: typ.ClassVar = {
 		'XYZ': {
-			'X': sockets.MaxwellBoundCondSocketDef(),
-			'Y': sockets.MaxwellBoundCondSocketDef(),
-			'Z': sockets.MaxwellBoundCondSocketDef(),
+			'X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
 		},
 		'±X | YZ': {
-			'+X': sockets.MaxwellBoundCondSocketDef(),
-			'-X': sockets.MaxwellBoundCondSocketDef(),
-			'Y': sockets.MaxwellBoundCondSocketDef(),
-			'Z': sockets.MaxwellBoundCondSocketDef(),
+			'+X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'-X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
 		},
 		'X | ±Y | Z': {
-			'X': sockets.MaxwellBoundCondSocketDef(),
-			'+Y': sockets.MaxwellBoundCondSocketDef(),
-			'-Y': sockets.MaxwellBoundCondSocketDef(),
-			'Z': sockets.MaxwellBoundCondSocketDef(),
+			'X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'+Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'-Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
 		},
 		'XY | ±Z': {
-			'X': sockets.MaxwellBoundCondSocketDef(),
-			'Y': sockets.MaxwellBoundCondSocketDef(),
-			'+Z': sockets.MaxwellBoundCondSocketDef(),
-			'-Z': sockets.MaxwellBoundCondSocketDef(),
+			'X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'+Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
+			'-Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
 		},
 		'±XY | Z': {
-			'+X': sockets.MaxwellBoundCondSocketDef(),
-			'-X': sockets.MaxwellBoundCondSocketDef(),
-			'+Y': sockets.MaxwellBoundCondSocketDef(),
-			'-Y': sockets.MaxwellBoundCondSocketDef(),
-			'Z': sockets.MaxwellBoundCondSocketDef(),
+			'+X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'-X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'+Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'-Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
 		},
 		'X | ±YZ': {
-			'X': sockets.MaxwellBoundCondSocketDef(),
-			'+Y': sockets.MaxwellBoundCondSocketDef(),
-			'-Y': sockets.MaxwellBoundCondSocketDef(),
-			'+Z': sockets.MaxwellBoundCondSocketDef(),
-			'-Z': sockets.MaxwellBoundCondSocketDef(),
+			'X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'+Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'-Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'+Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
+			'-Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
 		},
 		'±XYZ': {
-			'+X': sockets.MaxwellBoundCondSocketDef(),
-			'-X': sockets.MaxwellBoundCondSocketDef(),
-			'+Y': sockets.MaxwellBoundCondSocketDef(),
-			'-Y': sockets.MaxwellBoundCondSocketDef(),
-			'+Z': sockets.MaxwellBoundCondSocketDef(),
-			'-Z': sockets.MaxwellBoundCondSocketDef(),
+			'+X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'-X': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.X}),
+			'+Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'-Y': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Y}),
+			'+Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
+			'-Z': sockets.MaxwellBoundCondSocketDef(allow_axes={SSA.Z}),
 		},
 	}
 	output_sockets: typ.ClassVar = {

src/blender_maxwell/node_trees/maxwell_sim_nodes/nodes/simulations/sim_domain.py

@@ -35,27 +35,8 @@ class SimDomainNode(base.MaxwellSimNode):
 	}
 
 	####################
-	# - Event Methods
+	# - Events
 	####################
-	@events.computes_output_socket(
-		'Domain',
-		input_sockets={'Duration', 'Center', 'Size', 'Grid', 'Ambient Medium'},
-		unit_systems={'Tidy3DUnits': ct.UNITS_TIDY3D},
-		scale_input_sockets={
-			'Duration': 'Tidy3DUnits',
-			'Center': 'Tidy3DUnits',
-			'Size': 'Tidy3DUnits',
-		},
-	)
-	def compute_domain(self, input_sockets: dict, unit_systems) -> sp.Expr:
-		return {
-			'run_time': input_sockets['Duration'],
-			'center': input_sockets['Center'],
-			'size': input_sockets['Size'],
-			'grid_spec': input_sockets['Grid'],
-			'medium': input_sockets['Ambient Medium'],
-		}
-
 	@events.on_value_changed(
 		socket_name={'Center', 'Size'},
 		prop_name='preview_active',
@@ -91,6 +72,28 @@ class SimDomainNode(base.MaxwellSimNode):
 		if props['preview_active']:
 			managed_objs['mesh'].show_preview()
 
+	####################
+	# - Outputs
+	####################
+	@events.computes_output_socket(
+		'Domain',
+		input_sockets={'Duration', 'Center', 'Size', 'Grid', 'Ambient Medium'},
+		unit_systems={'Tidy3DUnits': ct.UNITS_TIDY3D},
+		scale_input_sockets={
+			'Duration': 'Tidy3DUnits',
+			'Center': 'Tidy3DUnits',
+			'Size': 'Tidy3DUnits',
+		},
+	)
+	def compute_domain(self, input_sockets, unit_systems) -> sp.Expr:
+		return {
+			'run_time': input_sockets['Duration'],
+			'center': input_sockets['Center'],
+			'size': input_sockets['Size'],
+			'grid_spec': input_sockets['Grid'],
+			'medium': input_sockets['Ambient Medium'],
+		}
+
 
 ####################
 # - Blender Registration

src/blender_maxwell/node_trees/maxwell_sim_nodes/sockets/base.py

@@ -373,6 +373,32 @@ class MaxwellSimSocket(bpy.types.NodeSocket):
 		"""
 		self.trigger_event(ct.FlowEvent.LinkChanged)
 
+	def remove_invalidated_links(self) -> None:
+		"""Reevaluates the capabilities of all socket links, and removes any that no longer match.
+
+		Links are removed with a simple `node_tree.links.remove()`, which directly emulates a user trying to remove the node link.
+		**Note** that all of the usual consent-semantics apply just the same as if the user had manually tried to remove the link.
+
+		Notes:
+			Called by nodes directly on their sockets, after altering any property that might influence the capabilities of that socket.
+
+			This prevents invalid use when the user alters a property, which **would** disallow adding a _new_ link identical to one that already exists.
+			In such a case, the existing (non-capability-respecting) link should be removed, as it has become invalid.
+		"""
+		node_tree = self.id_data
+		for link in self.links:
+			if not link.from_socket.capabilities.is_compatible_with(
+				link.to_socket.capabilities
+			):
+				log.error(
+					'Deleted link between "%s" (%s) and "%s" (%s) due to invalidated capabilities',
+					link.from_socket.bl_label,
+					link.from_socket.capabilities,
+					link.to_socket.bl_label,
+					link.to_socket.capabilities,
+				)
+				node_tree.links.remove(link)
+
 	####################
 	# - Event Chain
 	####################

src/blender_maxwell/node_trees/maxwell_sim_nodes/sockets/maxwell/bound_cond.py

@@ -26,6 +26,16 @@ class MaxwellBoundCondBLSocket(base.MaxwellSimSocket):
 	####################
 	default: ct.BoundCondType = bl_cache.BLField(ct.BoundCondType.Pml, prop_ui=True)
 
+	# Capabilities
+	## Allow a boundary condition compatible with any of the following axes.
+	allow_axes: set[ct.SimSpaceAxis] = bl_cache.BLField(
+		{ct.SimSpaceAxis.X, ct.SimSpaceAxis.Y, ct.SimSpaceAxis.Z},
+	)
+	## Present a boundary condition compatible with any of the following axes.
+	present_axes: set[ct.SimSpaceAxis] = bl_cache.BLField(
+		{ct.SimSpaceAxis.X, ct.SimSpaceAxis.Y, ct.SimSpaceAxis.Z},
+	)
+
 	####################
 	# - UI
 	####################
@@ -33,8 +43,17 @@ class MaxwellBoundCondBLSocket(base.MaxwellSimSocket):
 		col.prop(self, self.blfields['default'], text='')
 
 	####################
-	# - Computation of Default Value
+	# - FlowKind
 	####################
+	@property
+	def capabilities(self) -> ct.CapabilitiesFlow:
+		return ct.CapabilitiesFlow(
+			socket_type=self.socket_type,
+			active_kind=self.active_kind,
+			allow_any=self.allow_axes,
+			present_any=self.present_axes,
+		)
+
 	@property
 	def value(self) -> td.BoundaryEdge:
 		return self.default.tidy3d_boundary_edge
@@ -51,10 +70,23 @@ class MaxwellBoundCondSocketDef(base.SocketDef):
 	socket_type: ct.SocketType = ct.SocketType.MaxwellBoundCond
 
 	default: ct.BoundCondType = ct.BoundCondType.Pml
+	allow_axes: set[ct.SimSpaceAxis] = {
+		ct.SimSpaceAxis.X,
+		ct.SimSpaceAxis.Y,
+		ct.SimSpaceAxis.Z,
+	}
+	present_axes: set[ct.SimSpaceAxis] = {
+		ct.SimSpaceAxis.X,
+		ct.SimSpaceAxis.Y,
+		ct.SimSpaceAxis.Z,
+	}
 
 	def init(self, bl_socket: MaxwellBoundCondBLSocket) -> None:
 		bl_socket.default = self.default
 
+		bl_socket.allow_axes = self.allow_axes
+		bl_socket.present_axes = self.present_axes
+
 
 ####################
 # - Blender Registration

src/blender_maxwell/node_trees/maxwell_sim_nodes/sockets/maxwell/source.py

@@ -1,4 +1,3 @@
-
 from ... import contracts as ct
 from .. import base
 
@@ -16,6 +15,8 @@ class MaxwellSourceSocketDef(base.SocketDef):
 
 	is_list: bool = False
 
+	## TODO: capabilities() to require source sockets
+
 	def init(self, bl_socket: MaxwellSourceBLSocket) -> None:
 		if self.is_list:
 			bl_socket.active_kind = ct.FlowKind.Array