Geometry & Assembly

Base Geometry

BaseGeometry is the abstract base class for all geometric primitives and imported CAD models. It handles meshing, port definition, boundary detection, and persistence.

cavsim3d.geometry.base.BaseGeometry

Bases: ABC, TaggableMixin

Abstract base class for all geometries with tagging support.

Attributes

tag : ComponentTag Unique identifier for this component configuration compute_method : ComputeMethod Method for computing solution (NUMERIC, ANALYTICAL, SEMI_ANALYTICAL) custom_tag : str User-defined tag for identification

Source code in cavsim3d/geometry/base.py

class BaseGeometry(ABC, TaggableMixin):
    """
    Abstract base class for all geometries with tagging support.

    Attributes
    ----------
    tag : ComponentTag
        Unique identifier for this component configuration
    compute_method : ComputeMethod
        Method for computing solution (NUMERIC, ANALYTICAL, SEMI_ANALYTICAL)
    custom_tag : str
        User-defined tag for identification
    """

    _AXIS_MAP = {'X': X, 'Y': Y, 'Z': Z}
    _AXIS_ORDER = ['X', 'Y', 'Z']

    def __init__(self):
        self.geo = None
        self.mesh = None
        self.bc = None
        self._bc_explicitly_set = False
        self._ports = None
        self._boundaries = None

        # Tagging (from TaggableMixin)
        self._tag = None
        self._custom_tag = None
        self._compute_method = ComputeMethod.NUMERIC

        # History system (unified across all geometry types)
        self._history: List[dict] = []
        self._source_link: Optional[str] = None  # original file path
        self._source_hash: Optional[str] = None  # SHA-256 of source at save time

    def get_extreme_faces(self, axis: str = 'Z') -> Dict[str, Tuple[float, float, Any]]:
        """
        Identify the most extreme planar faces perpendicular to the axis.

        Returns
        -------
        dict
            'min': (position, area, face_obj)
            'max': (position, area, face_obj)
        """
        if self.geo is None:
            return {}

        try:
            axis_idx = self._AXIS_MAP[axis.upper()].index
        except (AttributeError, KeyError):
            axis_idx = {'X': 0, 'Y': 1, 'Z': 2}.get(axis.upper(), 2)

        extreme_min = {'pos': float('inf'), 'area': 0.0, 'face': None}
        extreme_max = {'pos': float('-inf'), 'area': 0.0, 'face': None}

        tolerance = 1e-4

        # First pass: find max area at each end
        try:
            for face in self.geo.faces:
                bb = face.bounding_box
                extent = bb[1][axis_idx] - bb[0][axis_idx]

                if extent < tolerance:
                    pos = (bb[0][axis_idx] + bb[1][axis_idx]) / 2

                    # Calculate proxy area from bounding box
                    axes = [0, 1, 2]
                    axes.remove(axis_idx)
                    area = (bb[1][axes[0]] - bb[0][axes[0]]) * (bb[1][axes[1]] - bb[0][axes[1]])

                    if pos < extreme_min['pos'] - tolerance:
                        extreme_min = {'pos': pos, 'area': area, 'face': face}
                    elif abs(pos - extreme_min['pos']) < tolerance:
                        if area > extreme_min['area']:
                            extreme_min = {'pos': pos, 'area': area, 'face': face}

                    if pos > extreme_max['pos'] + tolerance:
                        extreme_max = {'pos': pos, 'area': area, 'face': face}
                    elif abs(pos - extreme_max['pos']) < tolerance:
                        if area > extreme_max['area']:
                            extreme_max = {'pos': pos, 'area': area, 'face': face}

            # Second pass: ensure we didn't pick a tiny ghost face that's more extreme
            # but much smaller than the actual port face.
            max_face_area = max(extreme_min['area'], extreme_max['area'])
            if max_face_area > 0:
                refined_min = {'pos': float('inf'), 'area': 0.0, 'face': None}
                refined_max = {'pos': float('-inf'), 'area': 0.0, 'face': None}

                for face in self.geo.faces:
                    bb = face.bounding_box
                    if (bb[1][axis_idx] - bb[0][axis_idx]) < tolerance:
                        axes = [0, 1, 2]
                        axes.remove(axis_idx)
                        area = (bb[1][axes[0]] - bb[0][axes[0]]) * (bb[1][axes[1]] - bb[0][axes[1]])

                        if area > max_face_area * 0.05: # At least 5% of largest face
                            pos = (bb[0][axis_idx] + bb[1][axis_idx]) / 2
                            if pos < refined_min['pos']:
                                refined_min = {'pos': pos, 'area': area, 'face': face}
                            if pos > refined_max['pos']:
                                refined_max = {'pos': pos, 'area': area, 'face': face}

                if refined_min['face'] is not None: extreme_min = refined_min
                if refined_max['face'] is not None: extreme_max = refined_max

        except Exception as e:
            print(f"  [Warning] get_extreme_faces Exception: {e}")
            pass

        result = {}
        if extreme_min['face']:
            result['min'] = (extreme_min['pos'], extreme_min['area'], extreme_min['face'])
        if extreme_max['face']:
            result['max'] = (extreme_max['pos'], extreme_max['area'], extreme_max['face'])

        return result

    def get_physical_bounds(self, axis: str = 'Z') -> Tuple[float, float]:
        """
        Calculate physical extremes using the actual faces.
        """
        extremes = self.get_extreme_faces(axis)
        if 'min' in extremes and 'max' in extremes:
            return (extremes['min'][0], extremes['max'][0])

        # Fallback to general bounding box
        if self.geo is None:
            return (0.0, 1.0)

        try:
            axis_idx = self._AXIS_MAP[axis.upper()].index
        except (AttributeError, KeyError):
            axis_idx = {'X': 0, 'Y': 1, 'Z': 2}.get(axis.upper(), 2)

        bb = self.geo.bounding_box
        return (bb[0][axis_idx], bb[1][axis_idx])

    @abstractmethod
    def build(self) -> None:
        """Build the geometry."""
        pass

    def generate_mesh(self, maxh=None, curve_order: int = 3) -> Mesh:
        """Generate mesh from cavsim3d.geometry. Automatically calls build() if needed."""
        if self.geo is None:
            self.build()

        if maxh:
            self.mesh = Mesh(OCCGeometry(self.geo).GenerateMesh(maxh=maxh))
        else:
            self.mesh = Mesh(OCCGeometry(self.geo).GenerateMesh())

        self.mesh.Curve(curve_order)
        self.invalidate_tag()  # Mesh changed

        self._record('generate_mesh', maxh=maxh, curve_order=curve_order)
        return self.mesh

    def show(
        self,
        what: Literal["geometry", "mesh", "geo"] = "geometry",
        **kwargs
    ) -> None:
        """Display the geometry or mesh."""
        what = what.lower()

        if what in ("geometry", "geo"):
            if self.geo is None:
                self.build()
            scene = Draw(self.geo, **kwargs)
        elif what == "mesh":
            if self.mesh is None:
                raise ValueError("Mesh not generated. Call generate_mesh() first.")
            # Use NGSolve's Draw for NGSolve Mesh objects
            scene = Draw(self.mesh, **kwargs)
        else:
            raise ValueError(f"Invalid option '{what}'.")

        _display_webgui_fallback(scene)

    @property
    def ports(self) -> List[str]:
        """Get list of port boundary names."""
        if self._ports is None:
            self._ports = [b for b in self.mesh.GetBoundaries() if "port" in b]
        return self._ports

    @property
    def boundaries(self) -> List[str]:
        """Get all boundary names."""
        if self._boundaries is None:
            self._boundaries = list(self.mesh.GetBoundaries())
        return self._boundaries

    @property
    def n_ports(self) -> int:
        """Number of ports."""
        return len(self.ports)

    # === Caching helpers ===

    def get_cached_solution(self) -> Optional[CachedSolution]:
        """Get cached solution if available."""
        return get_global_cache().get(self.tag)

    def has_cached_solution(self) -> bool:
        """Check if a cached solution exists."""
        return self.tag in get_global_cache()

    def cache_solution(self, solution_data: Any, **metadata) -> None:
        """Store a solution in the cache."""
        get_global_cache().store(
            self.tag,
            solution_data,
            self.compute_method,
            **metadata
        )

    # === Keep all other methods from the original base.py ===
    # (define_ports, validate_boundaries, get_boundary_normal, etc.)

    def define_ports(self, **kwargs) -> 'BaseGeometry':
        """Define port faces by axis position."""
        # [Same implementation as before]
        # ...
        self.invalidate_tag()  # Ports changed
        return self

    def get_boundary_normal(self, boundary_label: str) -> Optional[np.ndarray]:
        """Calculate outward normal vector for a planar boundary face."""
        nhat = specialcf.normal(3)
        integral_n = Integrate(
            nhat, self.mesh, BND,
            definedon=self.mesh.Boundaries(boundary_label)
        )
        face_area = Integrate(
            1, self.mesh, BND,
            definedon=self.mesh.Boundaries(boundary_label)
        )
        if abs(face_area) < 1e-12:
            return None
        normal = -np.array(integral_n) / face_area
        return np.round(normal, decimals=6)

    def get_point_on_boundary(self, boundary_label: str) -> Optional[Tuple[float, ...]]:
        """Get coordinates of a vertex on the specified boundary."""
        for bel in self.mesh.Elements(BND):
            if bel.mat == boundary_label:
                vertex_index = bel.vertices[0]
                return tuple(self.mesh[vertex_index].point)
        return None

    def save_step(self, filename: str) -> None:
        """Save geometry to STEP file."""
        self._export('step', filename)

    def save_brep(self, filename: str) -> None:
        """Save geometry to BREP file."""
        self._export('brep', filename)

    def _export(self, format: Literal['step', 'brep', 'stl'], filename: str) -> None:
        """Export geometry to file."""
        if self.geo is None:
            raise ValueError("Geometry not built.")

        path = Path(filename)
        config = {
            'step': ('.step', lambda f: self.geo.WriteStep(str(f))),
            'brep': ('.brep', lambda f: self.geo.WriteBrep(str(f))),
            'stl': ('.stl', lambda f: self.geo.WriteSTL(str(f)))
        }
        ext, writer = config[format]
        if path.suffix.lower() != ext:
            path = path.with_suffix(ext)
        path.parent.mkdir(parents=True, exist_ok=True)
        writer(path)
        print(f"Saved to: {path}")

    def print_info(self) -> None:
        """Print detailed geometry information. Automatically builds if needed."""
        if self.geo is None:
            try:
                self.build()
            except Exception:
                pass

        class_name = self.__class__.__name__

        print("\n" + "=" * 70)
        print(f"{class_name} Geometry Information")
        print("=" * 70)

        print(f"Geometry type:          {class_name}")
        print(f"Compute method:         {self.compute_method}")
        print(f"Supports analytical:    {self.supports_analytical}")
        print(f"Boundary condition:     {self.bc}")

        print(f"\nComponent Tag:")
        print(f"  Full:                 {self.tag}")
        print(f"  Geometry hash:        {self.tag.geometry_hash[:16]}...")
        if self.custom_tag:
            print(f"  Custom tag:           {self.custom_tag}")

        has_cached = self.has_cached_solution()
        print(f"\nCache status:           {'CACHED' if has_cached else 'NOT CACHED'}")

        has_mesh = self.mesh is not None
        print(f"\nMesh generated:         {has_mesh}")
        if has_mesh:
            print(f"  Vertices:             {self.mesh.nv:,}")
            print(f"  Elements:             {self.mesh.ne:,}")
            print(f"  Ports:                {self.ports}")

        print("=" * 70)

    # === History system ===

    def _record(self, op: str, **params) -> None:
        """Append an operation to the history log."""
        entry = {
            'op': op,
            'timestamp': datetime.now().isoformat(),
        }
        # Filter out None values for cleaner history
        entry.update({k: v for k, v in params.items() if v is not None})
        self._history.append(entry)

    def get_history(self) -> List[dict]:
        """Return the full operation history."""
        return list(self._history)

    # === Geometry persistence ===

    @staticmethod
    def _file_hash(path: Path) -> str:
        """Compute SHA-256 hash of a file."""
        h = hashlib.sha256()
        with open(path, 'rb') as f:
            for chunk in iter(lambda: f.read(8192), b''):
                h.update(chunk)
        return h.hexdigest()

    def save_geometry(self, project_path: Path) -> None:
        """
        Save geometry files, history, and source link to the project.

        Creates ``project_path/geometry/`` containing:
        - ``source_model.<ext>`` — copy of the original source file
        - ``history.json`` — operation log + source link metadata
        """
        geo_dir = Path(project_path) / 'geometry'
        geo_dir.mkdir(parents=True, exist_ok=True)

        source_hash = None
        source_filename = None

        # Copy original source file if we have one
        if self._source_link is not None:
            src = Path(self._source_link)
            if src.exists():
                source_filename = f'source_model{src.suffix}'
                dest = geo_dir / source_filename
                shutil.copy2(str(src), str(dest))
                source_hash = self._file_hash(dest)
                self._source_hash = source_hash

        # Build history metadata
        meta = {
            'type': self.__class__.__name__,
            'module': self.__class__.__module__,
            'source_link': str(self._source_link) if self._source_link else None,
            'source_filename': source_filename,
            'source_hash': source_hash,
            'history': self._history,
        }

        with open(geo_dir / 'history.json', 'w') as f:
            json.dump(meta, f, indent=2, default=str)

    @classmethod
    def load_geometry(cls, project_path: Path) -> 'BaseGeometry':
        """
        Load geometry from a saved project by replaying the operation history.

        Reads ``project_path/geometry/history.json``, dispatches to the
        correct subclass, and reconstructs the geometry.
        """
        geo_dir = Path(project_path) / 'geometry'
        history_file = geo_dir / 'history.json'

        if not history_file.exists():
            raise FileNotFoundError(
                f"No geometry history found at {history_file}"
            )

        with open(history_file, 'r') as f:
            meta = json.load(f)

        geo_type = meta['type']
        history = meta['history']
        source_link = meta.get('source_link')
        source_filename = meta.get('source_filename')
        source_hash = meta.get('source_hash')

        # Resolve the geometry file path for replay
        # Use the project-local copy stored in geometry/
        source_file = geo_dir / source_filename if source_filename else None

        # Try to import the module stored in history to register the subclass
        if 'module' in meta:
            try:
                import importlib
                importlib.import_module(meta['module'])
            except (ImportError, KeyError):
                pass

        # Dispatch to the correct subclass
        subclass = cls._get_subclass(geo_type)
        if subclass is None:
            raise ValueError(
                f"Unknown geometry type '{geo_type}'. "
                f"Cannot reconstruct geometry."
            )

        geo = subclass._rebuild_from_history(history, project_path, source_file)
        geo._source_link = source_link
        geo._source_hash = source_hash
        geo._history = history

        # Check if source has changed since save
        geo._check_source_link(project_path)

        return geo

    @classmethod
    def _get_subclass(cls, type_name: str) -> Optional[type]:
        """Find a BaseGeometry subclass by name (searches all subclasses)."""
        # Lazy import to ensure all standard subclasses are registered if not already
        try:
            from . import primitives
            from . import importers
            from . import assembly
        except (ImportError, ValueError):
            pass

        def _search(klass):
            # Check the class itself
            if klass.__name__ == type_name:
                return klass
            # Search subclasses
            for sub in klass.__subclasses__():
                if sub.__name__ == type_name:
                    return sub
                found = _search(sub)
                if found:
                    return found
            return None

        return _search(cls)

    @classmethod
    def _rebuild_from_history(
        cls,
        history: List[dict],
        project_path: Path,
        source_file: Optional[Path] = None,
    ) -> 'BaseGeometry':
        """
        Reconstruct a geometry by replaying its history.

        Subclasses override this to handle their specific operations.
        """
        raise NotImplementedError(
            f"{cls.__name__} does not support history-based reconstruction."
        )

    def _check_source_link(self, project_path: Path) -> None:
        """
        Compare the project's geometry copy against the linked source.

        Interactive behaviour:
        - Source missing  → prompt to break link or keep it
        - Source changed → prompt to update (copy new), keep local, or break link
        - link_broken flag → skip all checks
        """
        if self._source_link is None:
            return  # No link — standalone project

        geo_dir = Path(project_path) / 'geometry'
        history_file = geo_dir / 'history.json'

        # Check if link was already broken
        if history_file.exists():
            with open(history_file, 'r') as f:
                meta = json.load(f)
            if meta.get('link_broken', False):
                return  # Link was broken — skip all checks

        source_path = Path(self._source_link)

        if not source_path.exists():
            print(f"\nLinked source geometry not found: {source_path}")
            print("  The original CAD file has been moved or deleted.")
            print("  The project's saved local copy will be used.\n")
            print("  Options:")
            print("    [1] Break link (stop checking in the future)")
            print("    [2] Keep link (ask again next time)")

            try:
                from cavsim3d.utils.io_utils import get_user_confirmation
                choice = input("  Choice [1/2]: ").strip()
                if choice == '1':
                    self._source_link = None
                    self._source_hash = None
                    # Write link_broken flag
                    if history_file.exists():
                        with open(history_file, 'r') as f:
                            meta = json.load(f)
                        meta['link_broken'] = True
                        meta['source_link'] = None
                        with open(history_file, 'w') as f:
                            json.dump(meta, f, indent=2, default=str)
                    print("  Link broken. Using local copy only.")
                else:
                    print("  Link preserved. Will check again next time.")
            except Exception:
                pass  # Non-interactive environment — just use local copy
            return

        # Source exists — check for changes
        try:
            current_hash = self._file_hash(source_path)
            if self._source_hash is not None and current_hash != self._source_hash:
                print(f"\n⚠ Linked source geometry has been modified: {source_path}")
                print("  The original CAD file has changed since the project was saved.")
                print("  Options:")
                print("    [1] Update (replace local copy with new file — invalidates results)")
                print("    [2] Keep local (ignore changes, use saved copy)")
                print("    [3] Break link (stop checking in the future)")

                try:
                    choice = input("  Choice [1/2/3]: ").strip()
                    if choice == '1':
                        # Copy new source into project
                        source_filename = f'source_model{source_path.suffix}'
                        dest = geo_dir / source_filename
                        shutil.copy2(str(source_path), str(dest))
                        self._source_hash = self._file_hash(dest)
                        self._update_link_in_history(geo_dir)
                        # Invalidate results
                        self._delete_project_results(project_path)
                        print("  Updated to new geometry. Results invalidated.")
                    elif choice == '3':
                        self._source_link = None
                        self._source_hash = None
                        if history_file.exists():
                            with open(history_file, 'r') as f:
                                meta = json.load(f)
                            meta['link_broken'] = True
                            meta['source_link'] = None
                            with open(history_file, 'w') as f:
                                json.dump(meta, f, indent=2, default=str)
                        print("  Link broken. Using local copy only.")
                    else:
                        print("  Keeping local copy. Original changes ignored.")
                except Exception:
                    pass  # Non-interactive — just use local copy
        except Exception as e:
            warnings.warn(f"Could not verify source geometry hash: {e}")


    def _update_link_in_history(self, geo_dir: Path) -> None:
        """Update the source_link and source_hash in history.json."""
        history_file = geo_dir / 'history.json'
        if history_file.exists():
            with open(history_file, 'r') as f:
                meta = json.load(f)
            meta['source_link'] = str(self._source_link) if self._source_link else None
            meta['source_hash'] = self._source_hash
            with open(history_file, 'w') as f:
                json.dump(meta, f, indent=2, default=str)

    @staticmethod
    def _delete_project_results(project_path: Path) -> None:
        """Delete all simulation results from a project directory."""
        result_paths = [
            'matrices.h5', 'snapshots.h5', 'fds',
            'fom', 'foms', 'roms', 'port_modes', 'eigenmode'
        ]
        for name in result_paths:
            p = Path(project_path) / name
            if p.is_file():
                p.unlink()
            elif p.is_dir():
                shutil.rmtree(p)

Functions

__init__()

Source code in cavsim3d/geometry/base.py

def __init__(self):
    self.geo = None
    self.mesh = None
    self.bc = None
    self._bc_explicitly_set = False
    self._ports = None
    self._boundaries = None

    # Tagging (from TaggableMixin)
    self._tag = None
    self._custom_tag = None
    self._compute_method = ComputeMethod.NUMERIC

    # History system (unified across all geometry types)
    self._history: List[dict] = []
    self._source_link: Optional[str] = None  # original file path
    self._source_hash: Optional[str] = None  # SHA-256 of source at save time

build() abstractmethod

Build the geometry.

Source code in cavsim3d/geometry/base.py

@abstractmethod
def build(self) -> None:
    """Build the geometry."""
    pass

generate_mesh(maxh=None, curve_order=3)

Generate mesh from cavsim3d.geometry. Automatically calls build() if needed.

Source code in cavsim3d/geometry/base.py

def generate_mesh(self, maxh=None, curve_order: int = 3) -> Mesh:
    """Generate mesh from cavsim3d.geometry. Automatically calls build() if needed."""
    if self.geo is None:
        self.build()

    if maxh:
        self.mesh = Mesh(OCCGeometry(self.geo).GenerateMesh(maxh=maxh))
    else:
        self.mesh = Mesh(OCCGeometry(self.geo).GenerateMesh())

    self.mesh.Curve(curve_order)
    self.invalidate_tag()  # Mesh changed

    self._record('generate_mesh', maxh=maxh, curve_order=curve_order)
    return self.mesh

define_ports(**kwargs)

Define port faces by axis position.

Source code in cavsim3d/geometry/base.py

def define_ports(self, **kwargs) -> 'BaseGeometry':
    """Define port faces by axis position."""
    # [Same implementation as before]
    # ...
    self.invalidate_tag()  # Ports changed
    return self

show(what='geometry', **kwargs)

Display the geometry or mesh.

Source code in cavsim3d/geometry/base.py

def show(
    self,
    what: Literal["geometry", "mesh", "geo"] = "geometry",
    **kwargs
) -> None:
    """Display the geometry or mesh."""
    what = what.lower()

    if what in ("geometry", "geo"):
        if self.geo is None:
            self.build()
        scene = Draw(self.geo, **kwargs)
    elif what == "mesh":
        if self.mesh is None:
            raise ValueError("Mesh not generated. Call generate_mesh() first.")
        # Use NGSolve's Draw for NGSolve Mesh objects
        scene = Draw(self.mesh, **kwargs)
    else:
        raise ValueError(f"Invalid option '{what}'.")

    _display_webgui_fallback(scene)

get_extreme_faces(axis='Z')

Identify the most extreme planar faces perpendicular to the axis.

Returns

dict 'min': (position, area, face_obj) 'max': (position, area, face_obj)

Source code in cavsim3d/geometry/base.py

def get_extreme_faces(self, axis: str = 'Z') -> Dict[str, Tuple[float, float, Any]]:
    """
    Identify the most extreme planar faces perpendicular to the axis.

    Returns
    -------
    dict
        'min': (position, area, face_obj)
        'max': (position, area, face_obj)
    """
    if self.geo is None:
        return {}

    try:
        axis_idx = self._AXIS_MAP[axis.upper()].index
    except (AttributeError, KeyError):
        axis_idx = {'X': 0, 'Y': 1, 'Z': 2}.get(axis.upper(), 2)

    extreme_min = {'pos': float('inf'), 'area': 0.0, 'face': None}
    extreme_max = {'pos': float('-inf'), 'area': 0.0, 'face': None}

    tolerance = 1e-4

    # First pass: find max area at each end
    try:
        for face in self.geo.faces:
            bb = face.bounding_box
            extent = bb[1][axis_idx] - bb[0][axis_idx]

            if extent < tolerance:
                pos = (bb[0][axis_idx] + bb[1][axis_idx]) / 2

                # Calculate proxy area from bounding box
                axes = [0, 1, 2]
                axes.remove(axis_idx)
                area = (bb[1][axes[0]] - bb[0][axes[0]]) * (bb[1][axes[1]] - bb[0][axes[1]])

                if pos < extreme_min['pos'] - tolerance:
                    extreme_min = {'pos': pos, 'area': area, 'face': face}
                elif abs(pos - extreme_min['pos']) < tolerance:
                    if area > extreme_min['area']:
                        extreme_min = {'pos': pos, 'area': area, 'face': face}

                if pos > extreme_max['pos'] + tolerance:
                    extreme_max = {'pos': pos, 'area': area, 'face': face}
                elif abs(pos - extreme_max['pos']) < tolerance:
                    if area > extreme_max['area']:
                        extreme_max = {'pos': pos, 'area': area, 'face': face}

        # Second pass: ensure we didn't pick a tiny ghost face that's more extreme
        # but much smaller than the actual port face.
        max_face_area = max(extreme_min['area'], extreme_max['area'])
        if max_face_area > 0:
            refined_min = {'pos': float('inf'), 'area': 0.0, 'face': None}
            refined_max = {'pos': float('-inf'), 'area': 0.0, 'face': None}

            for face in self.geo.faces:
                bb = face.bounding_box
                if (bb[1][axis_idx] - bb[0][axis_idx]) < tolerance:
                    axes = [0, 1, 2]
                    axes.remove(axis_idx)
                    area = (bb[1][axes[0]] - bb[0][axes[0]]) * (bb[1][axes[1]] - bb[0][axes[1]])

                    if area > max_face_area * 0.05: # At least 5% of largest face
                        pos = (bb[0][axis_idx] + bb[1][axis_idx]) / 2
                        if pos < refined_min['pos']:
                            refined_min = {'pos': pos, 'area': area, 'face': face}
                        if pos > refined_max['pos']:
                            refined_max = {'pos': pos, 'area': area, 'face': face}

            if refined_min['face'] is not None: extreme_min = refined_min
            if refined_max['face'] is not None: extreme_max = refined_max

    except Exception as e:
        print(f"  [Warning] get_extreme_faces Exception: {e}")
        pass

    result = {}
    if extreme_min['face']:
        result['min'] = (extreme_min['pos'], extreme_min['area'], extreme_min['face'])
    if extreme_max['face']:
        result['max'] = (extreme_max['pos'], extreme_max['area'], extreme_max['face'])

    return result

get_physical_bounds(axis='Z')

Calculate physical extremes using the actual faces.

Source code in cavsim3d/geometry/base.py

def get_physical_bounds(self, axis: str = 'Z') -> Tuple[float, float]:
    """
    Calculate physical extremes using the actual faces.
    """
    extremes = self.get_extreme_faces(axis)
    if 'min' in extremes and 'max' in extremes:
        return (extremes['min'][0], extremes['max'][0])

    # Fallback to general bounding box
    if self.geo is None:
        return (0.0, 1.0)

    try:
        axis_idx = self._AXIS_MAP[axis.upper()].index
    except (AttributeError, KeyError):
        axis_idx = {'X': 0, 'Y': 1, 'Z': 2}.get(axis.upper(), 2)

    bb = self.geo.bounding_box
    return (bb[0][axis_idx], bb[1][axis_idx])

get_boundary_normal(boundary_label)

Calculate outward normal vector for a planar boundary face.

Source code in cavsim3d/geometry/base.py

def get_boundary_normal(self, boundary_label: str) -> Optional[np.ndarray]:
    """Calculate outward normal vector for a planar boundary face."""
    nhat = specialcf.normal(3)
    integral_n = Integrate(
        nhat, self.mesh, BND,
        definedon=self.mesh.Boundaries(boundary_label)
    )
    face_area = Integrate(
        1, self.mesh, BND,
        definedon=self.mesh.Boundaries(boundary_label)
    )
    if abs(face_area) < 1e-12:
        return None
    normal = -np.array(integral_n) / face_area
    return np.round(normal, decimals=6)

get_point_on_boundary(boundary_label)

Get coordinates of a vertex on the specified boundary.

Source code in cavsim3d/geometry/base.py

def get_point_on_boundary(self, boundary_label: str) -> Optional[Tuple[float, ...]]:
    """Get coordinates of a vertex on the specified boundary."""
    for bel in self.mesh.Elements(BND):
        if bel.mat == boundary_label:
            vertex_index = bel.vertices[0]
            return tuple(self.mesh[vertex_index].point)
    return None

save_step(filename)

Save geometry to STEP file.

Source code in cavsim3d/geometry/base.py

def save_step(self, filename: str) -> None:
    """Save geometry to STEP file."""
    self._export('step', filename)

save_brep(filename)

Save geometry to BREP file.

Source code in cavsim3d/geometry/base.py

def save_brep(self, filename: str) -> None:
    """Save geometry to BREP file."""
    self._export('brep', filename)

print_info()

Print detailed geometry information. Automatically builds if needed.

Source code in cavsim3d/geometry/base.py

def print_info(self) -> None:
    """Print detailed geometry information. Automatically builds if needed."""
    if self.geo is None:
        try:
            self.build()
        except Exception:
            pass

    class_name = self.__class__.__name__

    print("\n" + "=" * 70)
    print(f"{class_name} Geometry Information")
    print("=" * 70)

    print(f"Geometry type:          {class_name}")
    print(f"Compute method:         {self.compute_method}")
    print(f"Supports analytical:    {self.supports_analytical}")
    print(f"Boundary condition:     {self.bc}")

    print(f"\nComponent Tag:")
    print(f"  Full:                 {self.tag}")
    print(f"  Geometry hash:        {self.tag.geometry_hash[:16]}...")
    if self.custom_tag:
        print(f"  Custom tag:           {self.custom_tag}")

    has_cached = self.has_cached_solution()
    print(f"\nCache status:           {'CACHED' if has_cached else 'NOT CACHED'}")

    has_mesh = self.mesh is not None
    print(f"\nMesh generated:         {has_mesh}")
    if has_mesh:
        print(f"  Vertices:             {self.mesh.nv:,}")
        print(f"  Elements:             {self.mesh.ne:,}")
        print(f"  Ports:                {self.ports}")

    print("=" * 70)

get_history()

Return the full operation history.

Source code in cavsim3d/geometry/base.py

def get_history(self) -> List[dict]:
    """Return the full operation history."""
    return list(self._history)

save_geometry(project_path)

Save geometry files, history, and source link to the project.

Creates project_path/geometry/ containing: - source_model.<ext> — copy of the original source file - history.json — operation log + source link metadata

Source code in cavsim3d/geometry/base.py

def save_geometry(self, project_path: Path) -> None:
    """
    Save geometry files, history, and source link to the project.

    Creates ``project_path/geometry/`` containing:
    - ``source_model.<ext>`` — copy of the original source file
    - ``history.json`` — operation log + source link metadata
    """
    geo_dir = Path(project_path) / 'geometry'
    geo_dir.mkdir(parents=True, exist_ok=True)

    source_hash = None
    source_filename = None

    # Copy original source file if we have one
    if self._source_link is not None:
        src = Path(self._source_link)
        if src.exists():
            source_filename = f'source_model{src.suffix}'
            dest = geo_dir / source_filename
            shutil.copy2(str(src), str(dest))
            source_hash = self._file_hash(dest)
            self._source_hash = source_hash

    # Build history metadata
    meta = {
        'type': self.__class__.__name__,
        'module': self.__class__.__module__,
        'source_link': str(self._source_link) if self._source_link else None,
        'source_filename': source_filename,
        'source_hash': source_hash,
        'history': self._history,
    }

    with open(geo_dir / 'history.json', 'w') as f:
        json.dump(meta, f, indent=2, default=str)

load_geometry(project_path) classmethod

Load geometry from a saved project by replaying the operation history.

Reads project_path/geometry/history.json, dispatches to the correct subclass, and reconstructs the geometry.

Source code in cavsim3d/geometry/base.py

@classmethod
def load_geometry(cls, project_path: Path) -> 'BaseGeometry':
    """
    Load geometry from a saved project by replaying the operation history.

    Reads ``project_path/geometry/history.json``, dispatches to the
    correct subclass, and reconstructs the geometry.
    """
    geo_dir = Path(project_path) / 'geometry'
    history_file = geo_dir / 'history.json'

    if not history_file.exists():
        raise FileNotFoundError(
            f"No geometry history found at {history_file}"
        )

    with open(history_file, 'r') as f:
        meta = json.load(f)

    geo_type = meta['type']
    history = meta['history']
    source_link = meta.get('source_link')
    source_filename = meta.get('source_filename')
    source_hash = meta.get('source_hash')

    # Resolve the geometry file path for replay
    # Use the project-local copy stored in geometry/
    source_file = geo_dir / source_filename if source_filename else None

    # Try to import the module stored in history to register the subclass
    if 'module' in meta:
        try:
            import importlib
            importlib.import_module(meta['module'])
        except (ImportError, KeyError):
            pass

    # Dispatch to the correct subclass
    subclass = cls._get_subclass(geo_type)
    if subclass is None:
        raise ValueError(
            f"Unknown geometry type '{geo_type}'. "
            f"Cannot reconstruct geometry."
        )

    geo = subclass._rebuild_from_history(history, project_path, source_file)
    geo._source_link = source_link
    geo._source_hash = source_hash
    geo._history = history

    # Check if source has changed since save
    geo._check_source_link(project_path)

    return geo

Properties

Bases: ABC, TaggableMixin

Abstract base class for all geometries with tagging support.

Attributes

tag : ComponentTag Unique identifier for this component configuration compute_method : ComputeMethod Method for computing solution (NUMERIC, ANALYTICAL, SEMI_ANALYTICAL) custom_tag : str User-defined tag for identification

Source code in cavsim3d/geometry/base.py

class BaseGeometry(ABC, TaggableMixin):
    """
    Abstract base class for all geometries with tagging support.

    Attributes
    ----------
    tag : ComponentTag
        Unique identifier for this component configuration
    compute_method : ComputeMethod
        Method for computing solution (NUMERIC, ANALYTICAL, SEMI_ANALYTICAL)
    custom_tag : str
        User-defined tag for identification
    """

    _AXIS_MAP = {'X': X, 'Y': Y, 'Z': Z}
    _AXIS_ORDER = ['X', 'Y', 'Z']

    def __init__(self):
        self.geo = None
        self.mesh = None
        self.bc = None
        self._bc_explicitly_set = False
        self._ports = None
        self._boundaries = None

        # Tagging (from TaggableMixin)
        self._tag = None
        self._custom_tag = None
        self._compute_method = ComputeMethod.NUMERIC

        # History system (unified across all geometry types)
        self._history: List[dict] = []
        self._source_link: Optional[str] = None  # original file path
        self._source_hash: Optional[str] = None  # SHA-256 of source at save time

    def get_extreme_faces(self, axis: str = 'Z') -> Dict[str, Tuple[float, float, Any]]:
        """
        Identify the most extreme planar faces perpendicular to the axis.

        Returns
        -------
        dict
            'min': (position, area, face_obj)
            'max': (position, area, face_obj)
        """
        if self.geo is None:
            return {}

        try:
            axis_idx = self._AXIS_MAP[axis.upper()].index
        except (AttributeError, KeyError):
            axis_idx = {'X': 0, 'Y': 1, 'Z': 2}.get(axis.upper(), 2)

        extreme_min = {'pos': float('inf'), 'area': 0.0, 'face': None}
        extreme_max = {'pos': float('-inf'), 'area': 0.0, 'face': None}

        tolerance = 1e-4

        # First pass: find max area at each end
        try:
            for face in self.geo.faces:
                bb = face.bounding_box
                extent = bb[1][axis_idx] - bb[0][axis_idx]

                if extent < tolerance:
                    pos = (bb[0][axis_idx] + bb[1][axis_idx]) / 2

                    # Calculate proxy area from bounding box
                    axes = [0, 1, 2]
                    axes.remove(axis_idx)
                    area = (bb[1][axes[0]] - bb[0][axes[0]]) * (bb[1][axes[1]] - bb[0][axes[1]])

                    if pos < extreme_min['pos'] - tolerance:
                        extreme_min = {'pos': pos, 'area': area, 'face': face}
                    elif abs(pos - extreme_min['pos']) < tolerance:
                        if area > extreme_min['area']:
                            extreme_min = {'pos': pos, 'area': area, 'face': face}

                    if pos > extreme_max['pos'] + tolerance:
                        extreme_max = {'pos': pos, 'area': area, 'face': face}
                    elif abs(pos - extreme_max['pos']) < tolerance:
                        if area > extreme_max['area']:
                            extreme_max = {'pos': pos, 'area': area, 'face': face}

            # Second pass: ensure we didn't pick a tiny ghost face that's more extreme
            # but much smaller than the actual port face.
            max_face_area = max(extreme_min['area'], extreme_max['area'])
            if max_face_area > 0:
                refined_min = {'pos': float('inf'), 'area': 0.0, 'face': None}
                refined_max = {'pos': float('-inf'), 'area': 0.0, 'face': None}

                for face in self.geo.faces:
                    bb = face.bounding_box
                    if (bb[1][axis_idx] - bb[0][axis_idx]) < tolerance:
                        axes = [0, 1, 2]
                        axes.remove(axis_idx)
                        area = (bb[1][axes[0]] - bb[0][axes[0]]) * (bb[1][axes[1]] - bb[0][axes[1]])

                        if area > max_face_area * 0.05: # At least 5% of largest face
                            pos = (bb[0][axis_idx] + bb[1][axis_idx]) / 2
                            if pos < refined_min['pos']:
                                refined_min = {'pos': pos, 'area': area, 'face': face}
                            if pos > refined_max['pos']:
                                refined_max = {'pos': pos, 'area': area, 'face': face}

                if refined_min['face'] is not None: extreme_min = refined_min
                if refined_max['face'] is not None: extreme_max = refined_max

        except Exception as e:
            print(f"  [Warning] get_extreme_faces Exception: {e}")
            pass

        result = {}
        if extreme_min['face']:
            result['min'] = (extreme_min['pos'], extreme_min['area'], extreme_min['face'])
        if extreme_max['face']:
            result['max'] = (extreme_max['pos'], extreme_max['area'], extreme_max['face'])

        return result

    def get_physical_bounds(self, axis: str = 'Z') -> Tuple[float, float]:
        """
        Calculate physical extremes using the actual faces.
        """
        extremes = self.get_extreme_faces(axis)
        if 'min' in extremes and 'max' in extremes:
            return (extremes['min'][0], extremes['max'][0])

        # Fallback to general bounding box
        if self.geo is None:
            return (0.0, 1.0)

        try:
            axis_idx = self._AXIS_MAP[axis.upper()].index
        except (AttributeError, KeyError):
            axis_idx = {'X': 0, 'Y': 1, 'Z': 2}.get(axis.upper(), 2)

        bb = self.geo.bounding_box
        return (bb[0][axis_idx], bb[1][axis_idx])

    @abstractmethod
    def build(self) -> None:
        """Build the geometry."""
        pass

    def generate_mesh(self, maxh=None, curve_order: int = 3) -> Mesh:
        """Generate mesh from cavsim3d.geometry. Automatically calls build() if needed."""
        if self.geo is None:
            self.build()

        if maxh:
            self.mesh = Mesh(OCCGeometry(self.geo).GenerateMesh(maxh=maxh))
        else:
            self.mesh = Mesh(OCCGeometry(self.geo).GenerateMesh())

        self.mesh.Curve(curve_order)
        self.invalidate_tag()  # Mesh changed

        self._record('generate_mesh', maxh=maxh, curve_order=curve_order)
        return self.mesh

    def show(
        self,
        what: Literal["geometry", "mesh", "geo"] = "geometry",
        **kwargs
    ) -> None:
        """Display the geometry or mesh."""
        what = what.lower()

        if what in ("geometry", "geo"):
            if self.geo is None:
                self.build()
            scene = Draw(self.geo, **kwargs)
        elif what == "mesh":
            if self.mesh is None:
                raise ValueError("Mesh not generated. Call generate_mesh() first.")
            # Use NGSolve's Draw for NGSolve Mesh objects
            scene = Draw(self.mesh, **kwargs)
        else:
            raise ValueError(f"Invalid option '{what}'.")

        _display_webgui_fallback(scene)

    @property
    def ports(self) -> List[str]:
        """Get list of port boundary names."""
        if self._ports is None:
            self._ports = [b for b in self.mesh.GetBoundaries() if "port" in b]
        return self._ports

    @property
    def boundaries(self) -> List[str]:
        """Get all boundary names."""
        if self._boundaries is None:
            self._boundaries = list(self.mesh.GetBoundaries())
        return self._boundaries

    @property
    def n_ports(self) -> int:
        """Number of ports."""
        return len(self.ports)

    # === Caching helpers ===

    def get_cached_solution(self) -> Optional[CachedSolution]:
        """Get cached solution if available."""
        return get_global_cache().get(self.tag)

    def has_cached_solution(self) -> bool:
        """Check if a cached solution exists."""
        return self.tag in get_global_cache()

    def cache_solution(self, solution_data: Any, **metadata) -> None:
        """Store a solution in the cache."""
        get_global_cache().store(
            self.tag,
            solution_data,
            self.compute_method,
            **metadata
        )

    # === Keep all other methods from the original base.py ===
    # (define_ports, validate_boundaries, get_boundary_normal, etc.)

    def define_ports(self, **kwargs) -> 'BaseGeometry':
        """Define port faces by axis position."""
        # [Same implementation as before]
        # ...
        self.invalidate_tag()  # Ports changed
        return self

    def get_boundary_normal(self, boundary_label: str) -> Optional[np.ndarray]:
        """Calculate outward normal vector for a planar boundary face."""
        nhat = specialcf.normal(3)
        integral_n = Integrate(
            nhat, self.mesh, BND,
            definedon=self.mesh.Boundaries(boundary_label)
        )
        face_area = Integrate(
            1, self.mesh, BND,
            definedon=self.mesh.Boundaries(boundary_label)
        )
        if abs(face_area) < 1e-12:
            return None
        normal = -np.array(integral_n) / face_area
        return np.round(normal, decimals=6)

    def get_point_on_boundary(self, boundary_label: str) -> Optional[Tuple[float, ...]]:
        """Get coordinates of a vertex on the specified boundary."""
        for bel in self.mesh.Elements(BND):
            if bel.mat == boundary_label:
                vertex_index = bel.vertices[0]
                return tuple(self.mesh[vertex_index].point)
        return None

    def save_step(self, filename: str) -> None:
        """Save geometry to STEP file."""
        self._export('step', filename)

    def save_brep(self, filename: str) -> None:
        """Save geometry to BREP file."""
        self._export('brep', filename)

    def _export(self, format: Literal['step', 'brep', 'stl'], filename: str) -> None:
        """Export geometry to file."""
        if self.geo is None:
            raise ValueError("Geometry not built.")

        path = Path(filename)
        config = {
            'step': ('.step', lambda f: self.geo.WriteStep(str(f))),
            'brep': ('.brep', lambda f: self.geo.WriteBrep(str(f))),
            'stl': ('.stl', lambda f: self.geo.WriteSTL(str(f)))
        }
        ext, writer = config[format]
        if path.suffix.lower() != ext:
            path = path.with_suffix(ext)
        path.parent.mkdir(parents=True, exist_ok=True)
        writer(path)
        print(f"Saved to: {path}")

    def print_info(self) -> None:
        """Print detailed geometry information. Automatically builds if needed."""
        if self.geo is None:
            try:
                self.build()
            except Exception:
                pass

        class_name = self.__class__.__name__

        print("\n" + "=" * 70)
        print(f"{class_name} Geometry Information")
        print("=" * 70)

        print(f"Geometry type:          {class_name}")
        print(f"Compute method:         {self.compute_method}")
        print(f"Supports analytical:    {self.supports_analytical}")
        print(f"Boundary condition:     {self.bc}")

        print(f"\nComponent Tag:")
        print(f"  Full:                 {self.tag}")
        print(f"  Geometry hash:        {self.tag.geometry_hash[:16]}...")
        if self.custom_tag:
            print(f"  Custom tag:           {self.custom_tag}")

        has_cached = self.has_cached_solution()
        print(f"\nCache status:           {'CACHED' if has_cached else 'NOT CACHED'}")

        has_mesh = self.mesh is not None
        print(f"\nMesh generated:         {has_mesh}")
        if has_mesh:
            print(f"  Vertices:             {self.mesh.nv:,}")
            print(f"  Elements:             {self.mesh.ne:,}")
            print(f"  Ports:                {self.ports}")

        print("=" * 70)

    # === History system ===

    def _record(self, op: str, **params) -> None:
        """Append an operation to the history log."""
        entry = {
            'op': op,
            'timestamp': datetime.now().isoformat(),
        }
        # Filter out None values for cleaner history
        entry.update({k: v for k, v in params.items() if v is not None})
        self._history.append(entry)

    def get_history(self) -> List[dict]:
        """Return the full operation history."""
        return list(self._history)

    # === Geometry persistence ===

    @staticmethod
    def _file_hash(path: Path) -> str:
        """Compute SHA-256 hash of a file."""
        h = hashlib.sha256()
        with open(path, 'rb') as f:
            for chunk in iter(lambda: f.read(8192), b''):
                h.update(chunk)
        return h.hexdigest()

    def save_geometry(self, project_path: Path) -> None:
        """
        Save geometry files, history, and source link to the project.

        Creates ``project_path/geometry/`` containing:
        - ``source_model.<ext>`` — copy of the original source file
        - ``history.json`` — operation log + source link metadata
        """
        geo_dir = Path(project_path) / 'geometry'
        geo_dir.mkdir(parents=True, exist_ok=True)

        source_hash = None
        source_filename = None

        # Copy original source file if we have one
        if self._source_link is not None:
            src = Path(self._source_link)
            if src.exists():
                source_filename = f'source_model{src.suffix}'
                dest = geo_dir / source_filename
                shutil.copy2(str(src), str(dest))
                source_hash = self._file_hash(dest)
                self._source_hash = source_hash

        # Build history metadata
        meta = {
            'type': self.__class__.__name__,
            'module': self.__class__.__module__,
            'source_link': str(self._source_link) if self._source_link else None,
            'source_filename': source_filename,
            'source_hash': source_hash,
            'history': self._history,
        }

        with open(geo_dir / 'history.json', 'w') as f:
            json.dump(meta, f, indent=2, default=str)

    @classmethod
    def load_geometry(cls, project_path: Path) -> 'BaseGeometry':
        """
        Load geometry from a saved project by replaying the operation history.

        Reads ``project_path/geometry/history.json``, dispatches to the
        correct subclass, and reconstructs the geometry.
        """
        geo_dir = Path(project_path) / 'geometry'
        history_file = geo_dir / 'history.json'

        if not history_file.exists():
            raise FileNotFoundError(
                f"No geometry history found at {history_file}"
            )

        with open(history_file, 'r') as f:
            meta = json.load(f)

        geo_type = meta['type']
        history = meta['history']
        source_link = meta.get('source_link')
        source_filename = meta.get('source_filename')
        source_hash = meta.get('source_hash')

        # Resolve the geometry file path for replay
        # Use the project-local copy stored in geometry/
        source_file = geo_dir / source_filename if source_filename else None

        # Try to import the module stored in history to register the subclass
        if 'module' in meta:
            try:
                import importlib
                importlib.import_module(meta['module'])
            except (ImportError, KeyError):
                pass

        # Dispatch to the correct subclass
        subclass = cls._get_subclass(geo_type)
        if subclass is None:
            raise ValueError(
                f"Unknown geometry type '{geo_type}'. "
                f"Cannot reconstruct geometry."
            )

        geo = subclass._rebuild_from_history(history, project_path, source_file)
        geo._source_link = source_link
        geo._source_hash = source_hash
        geo._history = history

        # Check if source has changed since save
        geo._check_source_link(project_path)

        return geo

    @classmethod
    def _get_subclass(cls, type_name: str) -> Optional[type]:
        """Find a BaseGeometry subclass by name (searches all subclasses)."""
        # Lazy import to ensure all standard subclasses are registered if not already
        try:
            from . import primitives
            from . import importers
            from . import assembly
        except (ImportError, ValueError):
            pass

        def _search(klass):
            # Check the class itself
            if klass.__name__ == type_name:
                return klass
            # Search subclasses
            for sub in klass.__subclasses__():
                if sub.__name__ == type_name:
                    return sub
                found = _search(sub)
                if found:
                    return found
            return None

        return _search(cls)

    @classmethod
    def _rebuild_from_history(
        cls,
        history: List[dict],
        project_path: Path,
        source_file: Optional[Path] = None,
    ) -> 'BaseGeometry':
        """
        Reconstruct a geometry by replaying its history.

        Subclasses override this to handle their specific operations.
        """
        raise NotImplementedError(
            f"{cls.__name__} does not support history-based reconstruction."
        )

    def _check_source_link(self, project_path: Path) -> None:
        """
        Compare the project's geometry copy against the linked source.

        Interactive behaviour:
        - Source missing  → prompt to break link or keep it
        - Source changed → prompt to update (copy new), keep local, or break link
        - link_broken flag → skip all checks
        """
        if self._source_link is None:
            return  # No link — standalone project

        geo_dir = Path(project_path) / 'geometry'
        history_file = geo_dir / 'history.json'

        # Check if link was already broken
        if history_file.exists():
            with open(history_file, 'r') as f:
                meta = json.load(f)
            if meta.get('link_broken', False):
                return  # Link was broken — skip all checks

        source_path = Path(self._source_link)

        if not source_path.exists():
            print(f"\nLinked source geometry not found: {source_path}")
            print("  The original CAD file has been moved or deleted.")
            print("  The project's saved local copy will be used.\n")
            print("  Options:")
            print("    [1] Break link (stop checking in the future)")
            print("    [2] Keep link (ask again next time)")

            try:
                from cavsim3d.utils.io_utils import get_user_confirmation
                choice = input("  Choice [1/2]: ").strip()
                if choice == '1':
                    self._source_link = None
                    self._source_hash = None
                    # Write link_broken flag
                    if history_file.exists():
                        with open(history_file, 'r') as f:
                            meta = json.load(f)
                        meta['link_broken'] = True
                        meta['source_link'] = None
                        with open(history_file, 'w') as f:
                            json.dump(meta, f, indent=2, default=str)
                    print("  Link broken. Using local copy only.")
                else:
                    print("  Link preserved. Will check again next time.")
            except Exception:
                pass  # Non-interactive environment — just use local copy
            return

        # Source exists — check for changes
        try:
            current_hash = self._file_hash(source_path)
            if self._source_hash is not None and current_hash != self._source_hash:
                print(f"\n⚠ Linked source geometry has been modified: {source_path}")
                print("  The original CAD file has changed since the project was saved.")
                print("  Options:")
                print("    [1] Update (replace local copy with new file — invalidates results)")
                print("    [2] Keep local (ignore changes, use saved copy)")
                print("    [3] Break link (stop checking in the future)")

                try:
                    choice = input("  Choice [1/2/3]: ").strip()
                    if choice == '1':
                        # Copy new source into project
                        source_filename = f'source_model{source_path.suffix}'
                        dest = geo_dir / source_filename
                        shutil.copy2(str(source_path), str(dest))
                        self._source_hash = self._file_hash(dest)
                        self._update_link_in_history(geo_dir)
                        # Invalidate results
                        self._delete_project_results(project_path)
                        print("  Updated to new geometry. Results invalidated.")
                    elif choice == '3':
                        self._source_link = None
                        self._source_hash = None
                        if history_file.exists():
                            with open(history_file, 'r') as f:
                                meta = json.load(f)
                            meta['link_broken'] = True
                            meta['source_link'] = None
                            with open(history_file, 'w') as f:
                                json.dump(meta, f, indent=2, default=str)
                        print("  Link broken. Using local copy only.")
                    else:
                        print("  Keeping local copy. Original changes ignored.")
                except Exception:
                    pass  # Non-interactive — just use local copy
        except Exception as e:
            warnings.warn(f"Could not verify source geometry hash: {e}")


    def _update_link_in_history(self, geo_dir: Path) -> None:
        """Update the source_link and source_hash in history.json."""
        history_file = geo_dir / 'history.json'
        if history_file.exists():
            with open(history_file, 'r') as f:
                meta = json.load(f)
            meta['source_link'] = str(self._source_link) if self._source_link else None
            meta['source_hash'] = self._source_hash
            with open(history_file, 'w') as f:
                json.dump(meta, f, indent=2, default=str)

    @staticmethod
    def _delete_project_results(project_path: Path) -> None:
        """Delete all simulation results from a project directory."""
        result_paths = [
            'matrices.h5', 'snapshots.h5', 'fds',
            'fom', 'foms', 'roms', 'port_modes', 'eigenmode'
        ]
        for name in result_paths:
            p = Path(project_path) / name
            if p.is_file():
                p.unlink()
            elif p.is_dir():
                shutil.rmtree(p)

ports property

Get list of port boundary names.

boundaries property

Get all boundary names.

n_ports property

Number of ports.

---

Assembly

Assembly manages multi-component geometries. Components are added sequentially and auto-aligned along a main axis. The assembly tracks connections (shared interfaces) between components and supports per-domain solving.

cavsim3d.geometry.assembly.Assembly

Bases: BaseGeometry

Assembly of geometry components with multi-solid output.

Creates a compound geometry where each component remains a separate solid, similar to split geometries. Supports automatic positioning, unified port naming, and tracks identical components for solver optimization.

Parameters

main_axis : str Primary axis for concatenation ('X', 'Y', or 'Z')

Examples

TESLA 9-cell cavity with identical midcells:

assembly = Assembly(main_axis='Z') assembly.add("endcell_l", endcell_l) for i in range(7): ... assembly.add("midcell", midcell, after="endcell_l" if i == 0 else "midcell") assembly.add("endcell_r", endcell_r, after="midcell") assembly.build()

Check what was created

assembly.print_info()

Shows: midcell_1, midcell_2, ..., midcell_7 (all identical)

Resulting structure: - Solids: endcell_l, midcell_1, midcell_2, ..., midcell_7, endcell_r - Ports: port1 (external), port2 (interface), ..., port10 (external)

Source code in cavsim3d/geometry/assembly.py

class Assembly(BaseGeometry):
    """
    Assembly of geometry components with multi-solid output.

    Creates a compound geometry where each component remains a separate solid,
    similar to split geometries. Supports automatic positioning, unified port
    naming, and tracks identical components for solver optimization.

    Parameters
    ----------
    main_axis : str
        Primary axis for concatenation ('X', 'Y', or 'Z')

    Examples
    --------
    **TESLA 9-cell cavity with identical midcells:**

    >>> assembly = Assembly(main_axis='Z')
    >>> assembly.add("endcell_l", endcell_l)
    >>> for i in range(7):
    ...     assembly.add("midcell", midcell, after="endcell_l" if i == 0 else "midcell")
    >>> assembly.add("endcell_r", endcell_r, after="midcell")
    >>> assembly.build()
    >>> 
    >>> # Check what was created
    >>> assembly.print_info()
    >>> # Shows: midcell_1, midcell_2, ..., midcell_7 (all identical)

    **Resulting structure:**
    - Solids: endcell_l, midcell_1, midcell_2, ..., midcell_7, endcell_r
    - Ports: port1 (external), port2 (interface), ..., port10 (external)
    """

    _AXIS_VEC = {'X': X, 'Y': Y, 'Z': Z}
    _AXIS_IDX = {'X': 0, 'Y': 1, 'Z': 2}

    def __init__(self, main_axis: Literal['X', 'Y', 'Z'] = 'Z'):
        super().__init__()
        self.main_axis = main_axis.upper()

        self._components: Dict[str, ComponentEntry] = {}
        self._component_order: List[str] = []
        self._connections: List[Connection] = []

        # Track base names for duplicate detection
        self._base_name_counts: Dict[str, int] = {}
        self._base_name_groups: Dict[str, List[str]] = {}  # base_name -> [key1, key2, ...]

        self._layout_computed = False
        self._is_built = False

        # Port and solid information (populated by build())
        self._port_info: Dict[str, Dict] = {}
        self._solid_info: Dict[str, Dict] = {}

        self._record('__init__', main_axis=main_axis)

    # =========================================================================
    # Component Management (with duplicate name support)
    # =========================================================================

    def add(
        self,
        name: str,
        geometry: Union[BaseGeometry, 'Assembly'],
        position: Optional[Tuple[float, float, float]] = None,
        rotation: Optional[Tuple[float, float, float]] = None,
        after: Optional[str] = None,
        before: Optional[str] = None,
        align_port: Optional[str] = None,
        **metadata
    ) -> 'Assembly':
        """
        Add a component (geometry or sub-assembly) to the assembly.

        Duplicate names are allowed and automatically suffixed with _1, _2, etc.
        Components with the same base name are tracked for solver optimization
        (compute once, reuse for identical geometries).

        Parameters
        ----------
        name : str
            Component name (duplicates allowed, will be auto-suffixed)
        geometry : BaseGeometry or Assembly
            Component to add (can be another Assembly!)
        position : tuple, optional
            Explicit position (x, y, z)
        rotation : tuple, optional
            Rotation angles in degrees (rx, ry, rz)
        after : str, optional
            Place after this component along main axis (can use base name)
        before : str, optional
            Place before this component along main axis (can use base name)
        align_port : str, optional
            Port to use for alignment
        **metadata
            Additional metadata

        Returns
        -------
        self : Assembly

        Examples
        --------
        >>> assembly.add("midcell", midcell)           # -> "midcell_1"
        >>> assembly.add("midcell", midcell)           # -> "midcell_2"  
        >>> assembly.add("midcell", midcell)           # -> "midcell_3"
        >>> assembly.get_identical_components()
        {'midcell': ['midcell_1', 'midcell_2', 'midcell_3']}
        """
        base_name = name

        # Generate unique key if name already exists
        if name in self._components or name in self._base_name_counts:
            # This is a duplicate name
            if base_name not in self._base_name_counts:
                # First duplicate - rename the original
                self._rename_first_instance(base_name)

            # Increment count and generate new key
            self._base_name_counts[base_name] += 1
            key = f"{base_name}_{self._base_name_counts[base_name]}"
        else:
            # First instance of this name
            key = name
            self._base_name_counts[base_name] = 0  # Will become 1 if duplicated

        # Track in base name groups
        if base_name not in self._base_name_groups:
            self._base_name_groups[base_name] = []
        self._base_name_groups[base_name].append(key)

        # Handle sub-assemblies
        if isinstance(geometry, Assembly):
            if not geometry._is_built:
                geometry.build()
        elif geometry.geo is None:
            # Auto-build if possible instead of raising error
            try:
                geometry.build()
            except Exception as e:
                raise ValueError(f"Geometry '{name}' not built and auto-build failed: {e}. Call build() first.")

        transform = Transform3D(
            translation=position or (0.0, 0.0, 0.0),
            rotation=rotation or (0.0, 0.0, 0.0)
        )

        entry = ComponentEntry(
            geometry=geometry,
            key=key,
            base_name=base_name,
            transform=transform,
            aligned_port=align_port,
            metadata=metadata
        )

        self._components[key] = entry

        # Handle ordering - resolve component and optional gap
        gap = 0.0
        if isinstance(after, tuple):
            ref_after, gap = after
        else:
            ref_after = after

        if isinstance(before, tuple):
            ref_before, gap = before
        else:
            ref_before = before

        resolved_after = self._resolve_component_ref(ref_after) if ref_after else None
        resolved_before = self._resolve_component_ref(ref_before) if ref_before else None

        if resolved_after is not None:
            idx = self._component_order.index(resolved_after) + 1
            self._component_order.insert(idx, key)

            self._connections.append(Connection(
                from_key=resolved_after,
                to_key=key,
                from_port='port2',
                to_port=align_port or 'port1',
                gap=gap
            ))

        elif resolved_before is not None:
            idx = self._component_order.index(resolved_before)
            self._component_order.insert(idx, key)

            self._connections.append(Connection(
                from_key=key,
                to_key=resolved_before,
                from_port=align_port or 'port2',
                to_port='port1',
                gap=-gap  # Before means moving current component backwards relative to target
            ))
        else:
            self._component_order.append(key)

        self._layout_computed = False
        self._is_built = False
        self.invalidate_tag()
        self._record(
            'add', 
            name=name, 
            geometry_type=type(geometry).__name__,
            geometry_history=geometry.get_history(),
            position=position,
            rotation=rotation,
            after=after,
            before=before,
            align_port=align_port,
            **metadata
        )

        return self

    def _rename_first_instance(self, base_name: str) -> None:
        """Rename the first instance of a component when duplicates are added."""
        if base_name not in self._components:
            return

        # The first instance needs to be renamed to base_name_1
        old_key = base_name
        new_key = f"{base_name}_1"

        # Update component entry
        entry = self._components[old_key]
        entry.key = new_key

        # Move in dict
        self._components[new_key] = entry
        del self._components[old_key]

        # Update order list
        idx = self._component_order.index(old_key)
        self._component_order[idx] = new_key

        # Update connections
        for conn in self._connections:
            if conn.from_key == old_key:
                conn.from_key = new_key
            if conn.to_key == old_key:
                conn.to_key = new_key

        # Update base name groups
        if base_name in self._base_name_groups:
            self._base_name_groups[base_name] = [
                new_key if k == old_key else k 
                for k in self._base_name_groups[base_name]
            ]

        # Set count to 1 (next will be _2)
        self._base_name_counts[base_name] = 1

    def _resolve_component_ref(self, ref: str) -> Optional[str]:
        """
        Resolve a component reference (base name or full key) to actual key.

        If ref is a base name with multiple instances, returns the last one.
        """
        if ref is None:
            return None

        # Direct match
        if ref in self._components:
            return ref

        # Try as base name - return last instance
        if ref in self._base_name_groups and self._base_name_groups[ref]:
            return self._base_name_groups[ref][-1]

        raise ValueError(f"Component '{ref}' not found")

    def get_identical_components(self) -> Dict[str, List[str]]:
        """
        Get groups of components that share the same base name (geometry).

        These components are candidates for solver optimization - compute
        the solution once and reuse for all identical instances.

        Returns
        -------
        dict
            {base_name: [key1, key2, ...]} for groups with more than one member
        """
        return {
            name: keys for name, keys in self._base_name_groups.items()
            if len(keys) > 1
        }

    def get_components_by_base_name(self, base_name: str) -> List[str]:
        """Get all component keys with a given base name."""
        return self._base_name_groups.get(base_name, [])

    # =========================================================================
    # Connection and Layout Methods (unchanged logic, updated for new structure)
    # =========================================================================

    def connect(
        self,
        name: str,
        geometry: Union[BaseGeometry, 'Assembly'],
        to_component: str,
        from_port: str = "port1",
        to_port: str = "port2",
        rotation: Optional[Tuple[float, float, float]] = None,
        gap: float = 0.0,
        **metadata
    ) -> 'Assembly':
        """Add component connected via ports to existing component."""
        resolved_to = self._resolve_component_ref(to_component)

        self.add(name, geometry, rotation=rotation, **metadata)

        # Get the key that was just added (might be suffixed)
        new_key = self._component_order[-1]

        self._record(
            'connect',
            name=name,
            geometry_type=type(geometry).__name__,
            geometry_history=geometry.get_history(),
            to_component=to_component,
            from_port=from_port,
            to_port=to_port,
            rotation=rotation,
            gap=gap,
            **metadata
        )

        return self

    def remove(self, ref: str) -> 'Assembly':
        """Remove a component by key or base name (removes all instances if base name)."""
        keys_to_remove = []

        if ref in self._components:
            keys_to_remove = [ref]
        elif ref in self._base_name_groups:
            keys_to_remove = list(self._base_name_groups[ref])
        else:
            raise KeyError(f"Component '{ref}' not found")

        for key in keys_to_remove:
            entry = self._components[key]
            base_name = entry.base_name

            del self._components[key]
            self._component_order.remove(key)
            self._connections = [
                c for c in self._connections 
                if c.from_key != key and c.to_key != key
            ]

            # Update base name tracking
            if base_name in self._base_name_groups:
                self._base_name_groups[base_name].remove(key)
                if not self._base_name_groups[base_name]:
                    del self._base_name_groups[base_name]
                    del self._base_name_counts[base_name]

        self._layout_computed = False
        self._is_built = False
        self.invalidate_tag()

        self._record('remove', ref=ref)

        return self

    def rotate(
        self, 
        ref: str, 
        angles: Tuple[float, float, float],
        center: Optional[Tuple[float, float, float]] = None
    ) -> 'Assembly':
        """Apply rotation to a component."""
        key = self._resolve_component_ref(ref)
        entry = self._components[key]
        current = entry.transform.rotation
        entry.transform = Transform3D(
            translation=entry.transform.translation,
            rotation=tuple(c + a for c, a in zip(current, angles)),
            rotation_center=center
        )
        self._layout_computed = False
        self._is_built = False
        return self

    def translate(self, ref: str, offset: Tuple[float, float, float]) -> 'Assembly':
        """Apply translation to a component."""
        key = self._resolve_component_ref(ref)
        entry = self._components[key]
        current = entry.transform.translation
        entry.transform = Transform3D(
            translation=tuple(c + o for c, o in zip(current, offset)),
            rotation=entry.transform.rotation,
            rotation_center=entry.transform.rotation_center
        )
        self._layout_computed = False
        self._is_built = False
        return self

    # =========================================================================
    # Index-based Access
    # =========================================================================

    def __getitem__(self, key: Union[str, int]) -> ComponentEntry:
        """Access component by key, base name, or index."""
        if isinstance(key, int):
            if key < 0 or key >= len(self._component_order):
                raise IndexError(f"Index {key} out of range")
            key = self._component_order[key]
        else:
            key = self._resolve_component_ref(key)

        return self._components[key]

    def __len__(self) -> int:
        return len(self._components)

    def __iter__(self):
        return iter(self._component_order)

    def __contains__(self, ref: str) -> bool:
        try:
            self._resolve_component_ref(ref)
            return True
        except (KeyError, ValueError):
            return False

    @property
    def keys(self) -> List[str]:
        return list(self._component_order)

    @property
    def components(self) -> Dict[str, ComponentEntry]:
        return dict(self._components)

    # =========================================================================
    # Layout Computation (unchanged)
    # =========================================================================

    def compute_layout(self) -> 'Assembly':
        """Compute positions based on connections."""
        if self._layout_computed:
            return self

        if not self._components:
            self._layout_computed = True
            return self

        positioned = set()
        first_key = self._component_order[0]
        positioned.add(first_key)

        for conn in self._connections:
            if conn.to_key in positioned and conn.from_key in positioned:
                continue
            if conn.from_key in positioned:
                self._position_connected(conn)
                positioned.add(conn.to_key)
            elif conn.to_key in positioned:
                # Reverse connection positioning not yet implemented for all cases
                pass

        # Safety for unpositioned components
        for key in self._component_order:
            if key not in positioned:
                if positioned:
                    last = list(positioned)[-1]
                    self._position_after(key, last)
                positioned.add(key)

        self._layout_computed = True
        return self

    def _get_port_position(
        self, 
        entry: ComponentEntry, 
        port_name: str
    ) -> Optional[Tuple[float, float, float]]:
        """Get port position in local coordinates using physical bounds."""
        # Refresh physical bounds for better precision if possible
        if not isinstance(entry.geometry, Assembly) and entry.geometry.geo is not None:
            axis = self.main_axis
            z_min, z_max = entry.geometry.get_physical_bounds('Z')
            x_min, x_max = entry.geometry.get_physical_bounds('X')
            y_min, y_max = entry.geometry.get_physical_bounds('Y')
            pmin = (x_min, y_min, z_min)
            pmax = (x_max, y_max, z_max)
        else:
            if entry.original_bounds is None:
                return None
            pmin, pmax = entry.original_bounds

        axis_idx = self._AXIS_IDX[self.main_axis]
        center = [(pmin[i] + pmax[i]) / 2 for i in range(3)]

        if 'port1' in port_name.lower() or 'min' in port_name.lower():
            center[axis_idx] = pmin[axis_idx]
        elif 'port2' in port_name.lower() or 'max' in port_name.lower():
            center[axis_idx] = pmax[axis_idx]

        return tuple(center)

    def _position_connected(self, conn: Connection) -> None:
        """Position component based on port connection."""
        from_entry = self._components[conn.from_key]
        to_entry = self._components[conn.to_key]

        from_port_pos = self._get_port_position(from_entry, conn.from_port)
        to_port_pos = self._get_port_position(to_entry, conn.to_port)

        if from_port_pos is None or to_port_pos is None:
            self._position_after(conn.to_key, conn.from_key)
            return

        axis_idx = self._AXIS_IDX[self.main_axis]

        from_world = tuple(
            from_port_pos[i] + from_entry.transform.translation[i]
            for i in range(3)
        )

        translation = [0.0, 0.0, 0.0]

        for i in range(3):
            if i == axis_idx:
                translation[i] = from_world[i] - to_port_pos[i] + conn.gap
            else:
                from_center = (from_entry.original_bounds[0][i] + 
                             from_entry.original_bounds[1][i]) / 2
                to_center = (to_entry.original_bounds[0][i] + 
                           to_entry.original_bounds[1][i]) / 2
                translation[i] = (from_center + from_entry.transform.translation[i] 
                                 - to_center)

        to_entry.transform = Transform3D(
            translation=tuple(translation),
            rotation=to_entry.transform.rotation,
            rotation_center=to_entry.transform.rotation_center
        )

    def _position_after(self, key: str, after_key: str, gap: float = 0.0) -> None:
        """Position component after another along main axis using bounding boxes."""
        entry = self._components[key]
        after_entry = self._components[after_key]

        axis_idx = self._AXIS_IDX[self.main_axis]

        # World max of previous component
        pmin_after, pmax_after = after_entry.original_bounds
        world_max_after = pmax_after[axis_idx] + after_entry.transform.translation[axis_idx]

        # Local min of current component
        pmin_self, pmax_self = entry.original_bounds
        local_min_self = pmin_self[axis_idx]

        translation = [0.0, 0.0, 0.0]

        for i in range(3):
            if i == axis_idx:
                translation[i] = world_max_after - local_min_self + gap
            else:
                after_center = (pmin_after[i] + pmax_after[i]) / 2 + after_entry.transform.translation[i]
                self_center = (pmin_self[i] + pmax_self[i]) / 2
                translation[i] = after_center - self_center

        entry.transform = Transform3D(
            translation=tuple(translation),
            rotation=entry.transform.rotation,
            rotation_center=entry.transform.rotation_center
        )

    def _position_before(self, key: str, before_key: str, gap: float = 0.0) -> None:
        """Position component before another along main axis using bounding boxes."""
        entry = self._components[key]
        before_entry = self._components[before_key]

        axis_idx = self._AXIS_IDX[self.main_axis]

        # World min of next component
        pmin_before, pmax_before = before_entry.original_bounds
        world_min_before = pmin_before[axis_idx] + before_entry.transform.translation[axis_idx]

        # Local max of current component
        pmin_self, pmax_self = entry.original_bounds
        local_max_self = pmax_self[axis_idx]

        translation = [0.0, 0.0, 0.0]

        for i in range(3):
            if i == axis_idx:
                translation[i] = world_min_before - local_max_self - gap
            else:
                before_center = (pmin_before[i] + pmax_before[i]) / 2 + before_entry.transform.translation[i]
                self_center = (pmin_self[i] + pmax_self[i]) / 2
                translation[i] = before_center - self_center

        entry.transform = Transform3D(
            translation=tuple(translation),
            rotation=entry.transform.rotation,
            rotation_center=entry.transform.rotation_center
        )

    # =========================================================================
    # Build - Multi-solid Assembly
    # =========================================================================

    def build(self) -> None:
        """
        Build multi-solid geometry from all components.

        Creates a compound where each component is a separate solid with:
        - Solid named after component key (e.g., 'midcell_1', 'midcell_2')
        - Sequential port naming (port1, port2, ...) along main axis
        - Interface ports (where solids meet) shared between adjacent components
        - External ports at assembly ends
        - All other faces named 'wall'
        """
        if not self._components:
            raise ValueError("No components in assembly")

        self.compute_layout()

        # NEW: Snap interfaces to compensate for CAD inaccuracies
        self._snap_interfaces()

        shapes = []
        # Build sub-assemblies first
        for entry in self._components.values():
            if entry.is_assembly and not entry.geometry._is_built:
                entry.geometry.build()

        # Collect transformed shapes
        shapes = []
        shape_keys = []  # Track which key each shape belongs to

        for key in self._component_order:
            entry = self._components[key]
            shape = entry.geometry.geo
            transformed = self._apply_transform(shape, entry.transform)
            shapes.append(transformed)
            shape_keys.append(key)

        # To prevent merging identical/co-planar solids, we assign unique 
        # temporary materials and use OCCGeometry(shapes).shape which is
        # often more robust than Glue() for preserving solid identities.
        if len(shapes) > 1:
            for i, shape in enumerate(shapes):
                shape.mat(f"__temp_solid_{i}")

        # Create multi-solid compound
        if len(shapes) == 1:
            self.geo = shapes[0]
        else:
            # OCCGeometry(list) performs gluing but preserves solid domains better
            occ_geo = OCCGeometry(shapes)
            self.geo = occ_geo.shape

        # Name solids and setup boundaries
        self._name_solids(shape_keys)
        self._setup_boundaries()

        self._is_built = True
        self._record('build')

    def generate_mesh(self, maxh=None, curve_order: int = 3) -> Mesh:
        """
        Generate mesh with support for per-component refinement.
        """
        if self.geo is None:
            self.build()

        # Apply per-component maxh if available
        solids = list(self.geo.solids)
        if len(solids) == len(self._component_order):
            for i, key in enumerate(self._component_order):
                entry = self._components[key]
                # Try to get maxh from component geometry
                comp_maxh = getattr(entry.geometry, 'maxh', None)
                if comp_maxh:
                    solids[i].maxh = comp_maxh

        # Call base generate_mesh which will create OCCGeometry(self.geo)
        # and respect the per-solid maxh settings.
        return super().generate_mesh(maxh=maxh, curve_order=curve_order)

    def _name_solids(self, shape_keys: List[str]) -> None:
        """Name each solid in the geometry based on component keys.

        Each assembly component may contain multiple mesh solids (e.g. after
        PEC subtraction).  All solids belonging to the same component are
        given the same material name so the solver treats each component as
        a single domain.
        """
        self._solid_info = {}
        axis_idx = self._AXIS_IDX[self.main_axis]

        try:
            solids = list(self.geo.solids)
        except AttributeError:
            # Single solid
            key = shape_keys[0]
            entry = self._components[key]
            self.geo.mat(key)
            self._solid_info[key] = {
                'material': key,
                'base_name': entry.base_name,
                'position': 0,
                'bounds': None,
                'is_identical': False,
            }
            return

        # Compute bounding box for each component along the main axis
        component_ranges = []
        for key in shape_keys:
            entry = self._components[key]
            trf = entry.transform
            offset = trf.translation[axis_idx]
            # Get the component's own extent along main axis
            geo = entry.geometry.geo
            try:
                comp_solids = list(geo.solids)
            except AttributeError:
                comp_solids = [geo]
            lo = min(s.bounding_box[0][axis_idx] for s in comp_solids) + offset
            hi = max(s.bounding_box[1][axis_idx] for s in comp_solids) + offset
            component_ranges.append((lo, hi, key))

        # Sort components along main axis
        component_ranges.sort(key=lambda x: (x[0] + x[1]) / 2)

        # Assign each mesh solid to the component whose range contains it.
        # Solid names are prefixed with the component key so the solver can
        # map mesh materials back to their owning component while preserving
        # per-solid material distinctions for CoefficientFunction assembly.
        def get_solid_center(solid):
            bb = solid.bounding_box
            return (bb[0][axis_idx] + bb[1][axis_idx]) / 2

        # domain_materials: component_key -> [mesh_material_name, ...]
        self._domain_materials: Dict[str, List[str]] = {key: [] for _, _, key in component_ranges}

        for idx, solid in enumerate(solids):
            center = get_solid_center(solid)
            best_key = None
            best_dist = float('inf')
            for lo, hi, key in component_ranges:
                if lo - 0.01 <= center <= hi + 0.01:
                    best_key = key
                    break
                mid = (lo + hi) / 2
                dist = abs(center - mid)
                if dist < best_dist:
                    best_dist = dist
                    best_key = key

            # Prefix original solid name with component key for uniqueness
            orig_name = solid.name if solid.name else f"solid"
            # Strip temporary assembly naming
            if orig_name.startswith('__temp_solid_'):
                orig_name = f"solid"
            mat_name = f"{best_key}/{orig_name}"
            solid.mat(mat_name)
            self._domain_materials[best_key].append(mat_name)

        # Print domain material mapping
        print(f"\nAssembly domain-material mapping:")
        for key, mats in self._domain_materials.items():
            print(f"  {key}: {mats}")

        # Store component-level info
        for lo, hi, key in component_ranges:
            entry = self._components[key]
            self._solid_info[key] = {
                'material': key,
                'base_name': entry.base_name,
                'position': (lo + hi) / 2,
                'bounds': (lo, hi),
                'is_identical': len(self._base_name_groups.get(entry.base_name, [])) > 1,
                'mesh_materials': self._domain_materials[key],
            }

    def _snap_interfaces(self, tolerance: float = 1e-1) -> None: # 0.1 meter = 10cm tolerance (very aggressive)
        """
        Perform a bit-perfect alignment pass to snap faces together.
        """
        if not self._component_order or len(self._component_order) < 2:
            return

        axis_idx = self._AXIS_IDX[self.main_axis]

        # We follow the component order and snap each to its predecessor
        for i in range(1, len(self._component_order)):
            to_key = self._component_order[i]
            from_key = self._component_order[i-1]

            from_entry = self._components[from_key]
            to_entry = self._components[to_key]

            if from_entry.geometry.geo is None or to_entry.geometry.geo is None:
                continue

            # Get BIT-PRECISE physical extremes from extreme faces
            from_extremes = from_entry.geometry.get_extreme_faces(self.main_axis)
            to_extremes = to_entry.geometry.get_extreme_faces(self.main_axis)

            if 'max' not in from_extremes or 'min' not in to_extremes:
                # Fallback to physical bounds if no extreme face identified
                _, pmax_from = from_entry.geometry.get_physical_bounds(self.main_axis)
                pmin_to, _ = to_entry.geometry.get_physical_bounds(self.main_axis)
            else:
                pmax_from = from_extremes['max'][0]
                pmin_to = to_extremes['min'][0]

            world_max_from = pmax_from + from_entry.transform.translation[axis_idx]
            world_min_to = pmin_to + to_entry.transform.translation[axis_idx]

            gap = world_min_to - world_max_from

            # SNAP SAFETY: If the gap is massive, something is wrong with face detection
            # We don't want to fly components across the world.
            # Usually a gap should be < 10% of component size for snapping.
            comp_size = max(from_entry.size[axis_idx], to_entry.size[axis_idx])
            if abs(gap) > comp_size * 0.5:
                # print(f"SNAP REJECTED: Gap too large ({gap:.3e}) for {to_key}")
                continue

            # Snap them if they are reasonably close
            if abs(gap) < tolerance:
                trans = list(to_entry.transform.translation)
                trans[axis_idx] -= gap
                to_entry.transform.translation = tuple(trans)
                self._layout_computed = True

    def _name_solids_by_position(self, solids: list) -> None:
        """Fallback: name solids by their position along main axis."""
        axis_idx = self._AXIS_IDX[self.main_axis]

        def get_solid_position(solid):
            bb = solid.bounding_box
            return (bb[0][axis_idx] + bb[1][axis_idx]) / 2

        solids_sorted = sorted(solids, key=get_solid_position)

        for i, solid in enumerate(solids_sorted):
            name = f"cell_{i + 1}"
            solid.mat(name)
            bb = solid.bounding_box
            self._solid_info[name] = {
                'material': name,
                'base_name': name,
                'position': get_solid_position(solid),
                'bounds': (tuple(bb[0]), tuple(bb[1])),
                'is_identical': False
            }

    def _setup_boundaries(self) -> None:
        """
        Setup boundary conditions with unified sequential port naming.

        - All flat faces perpendicular to main axis are ports
        - Faces at the same position share the same port name (interfaces)
        - Port numbers increase along the positive main axis direction
        - Non-port faces are named 'wall'
        """
        if self.geo is None:
            return

        axis_idx = self._AXIS_IDX[self.main_axis]
        tolerance = 1e-3 # 1mm grouping tolerance for ports

        # Step 1: Name all faces 'wall' first
        self._name_all_faces_wall()

        # Step 2: Identify all port faces and their positions
        port_faces = []

        try:
            for face in self.geo.faces:
                try:
                    fbb = face.bounding_box
                    face_extent = fbb[1][axis_idx] - fbb[0][axis_idx]

                    if face_extent < tolerance:
                        face_pos = (fbb[0][axis_idx] + fbb[1][axis_idx]) / 2
                        port_faces.append((face_pos, face))
                except Exception:
                    pass
        except Exception as e:
            print(f"Warning: Error identifying faces: {e}")
            self.bc = 'default'
            self._bc_explicitly_set = True
            return

        if not port_faces:
            print("Warning: No port faces found")
            self.bc = 'default'
            self._bc_explicitly_set = True
            return

        # Step 3: Group faces by position
        port_faces.sort(key=lambda x: x[0])

        port_groups = []
        for pos, face in port_faces:
            matched = False
            for group in port_groups:
                if abs(group['position'] - pos) < tolerance:
                    group['faces'].append(face)
                    matched = True
                    break

            if not matched:
                port_groups.append({
                    'position': pos,
                    'faces': [face]
                })

        # Step 4: Sort groups by position and assign sequential port names
        port_groups.sort(key=lambda g: g['position'])

        self._port_info = {}

        for port_num, group in enumerate(port_groups, start=1):
            port_name = f'port{port_num}'

            is_external = len(group['faces']) == 1
            is_interface = len(group['faces']) > 1

            for face in group['faces']:
                face.name = port_name
                face.col = (1, 0, 0)

            # Find which components connect at this port
            connected_components = self._find_components_at_position(
                group['position'], tolerance
            )

            self._port_info[port_name] = {
                'position': group['position'],
                'num_faces': len(group['faces']),
                'type': 'external' if is_external else 'interface',
                'axis': self.main_axis,
                'connected_components': connected_components
            }

        # Summary
        n_external = sum(1 for p in self._port_info.values() if p['type'] == 'external')
        n_interface = sum(1 for p in self._port_info.values() if p['type'] == 'interface')

        print(f"Port naming complete:")
        print(f"  Total ports: {len(self._port_info)}")
        print(f"  External ports: {n_external}")
        print(f"  Interface ports: {n_interface}")

        self.bc = 'default'
        self._bc_explicitly_set = True

    def _name_all_faces_wall(self) -> None:
        """Name every face 'wall' as default."""
        if self.geo is None:
            return

        try:
            for solid in self.geo.solids:
                for face in solid.faces:
                    face.name = 'default'
        except AttributeError:
            try:
                for face in self.geo.faces:
                    face.name = 'default'
            except AttributeError:
                pass

    def get_material(self, domain_name: str) -> dict:
        """Get material properties for a mesh material, delegating to its component.

        Mesh materials for assemblies are prefixed with the component key,
        e.g. ``cell1/ceramic_1``.  This method strips the prefix and queries
        the owning component's geometry for the material properties.
        """
        defaults = {"eps_r": 1.0, "mu_r": 1.0, "sigma": 0.0, "tan_delta": 0.0}

        # Try to find the owning component via the domain_materials mapping
        dm = getattr(self, '_domain_materials', {})
        for comp_key, mat_names in dm.items():
            if domain_name in mat_names:
                entry = self._components[comp_key]
                # Strip component prefix to get the original solid name
                orig_name = domain_name[len(comp_key) + 1:] if domain_name.startswith(comp_key + '/') else domain_name
                if hasattr(entry.geometry, 'get_material'):
                    return entry.geometry.get_material(orig_name)
                return dict(defaults)

        # Fallback: try component key as domain_name
        if domain_name in self._components:
            return dict(defaults)

        return dict(defaults)

    def _find_components_at_position(
        self, 
        position: float, 
        tolerance: float
    ) -> List[Tuple[str, str]]:
        """Find which components have ports at a given position."""
        axis_idx = self._AXIS_IDX[self.main_axis]
        connected = []

        for key in self._component_order:
            entry = self._components[key]
            if entry.original_bounds is None:
                continue

            pmin, pmax = entry.original_bounds
            t = entry.transform.translation

            world_min = pmin[axis_idx] + t[axis_idx]
            world_max = pmax[axis_idx] + t[axis_idx]

            if abs(world_min - position) < tolerance:
                connected.append((key, 'port1'))
            if abs(world_max - position) < tolerance:
                connected.append((key, 'port2'))

        return connected

    def _apply_transform(self, shape, transform: Transform3D):
        """Apply transformation to OCC shape."""
        if transform.is_identity():
            return shape

        result = shape

        rx, ry, rz = transform.rotation
        if abs(rx) > 1e-12 or abs(ry) > 1e-12 or abs(rz) > 1e-12:
            if transform.rotation_center:
                center = transform.rotation_center
            else:
                try:
                    bb = shape.bounding_box
                    center = tuple((bb[0][i] + bb[1][i]) / 2 for i in range(3))
                except:
                    center = (0, 0, 0)

            if abs(rz) > 1e-12:
                result = result.Rotate(Axes(center, Z), rz)
            if abs(ry) > 1e-12:
                result = result.Rotate(Axes(center, Y), ry)
            if abs(rx) > 1e-12:
                result = result.Rotate(Axes(center, X), rx)

        tx, ty, tz = transform.translation
        if abs(tx) > 1e-12 or abs(ty) > 1e-12 or abs(tz) > 1e-12:
            result = result.Move((tx, ty, tz))

        return result

    # =========================================================================
    # Port and Solid Information
    # =========================================================================

    def get_port_info(self) -> Dict[str, Dict]:
        """Get information about all ports in the assembly."""
        return dict(self._port_info) if self._port_info else {}

    def get_solid_info(self) -> Dict[str, Dict]:
        """Get information about all solids in the assembly."""
        return dict(self._solid_info) if self._solid_info else {}

    def get_external_ports(self) -> List[str]:
        """Get list of external port names."""
        return [name for name, info in self._port_info.items() 
                if info['type'] == 'external']

    def get_interface_ports(self) -> List[str]:
        """Get list of interface port names."""
        return [name for name, info in self._port_info.items() 
                if info['type'] == 'interface']

    def print_port_info(self) -> None:
        """Print detailed port information."""
        if not self._port_info:
            print("No port information available. Call build() first.")
            return

        print("\n" + "=" * 70)
        print("ASSEMBLY PORT MAP")
        print("=" * 70)
        print(f"Main axis: {self.main_axis}")
        print(f"Total ports: {len(self._port_info)}")
        print("-" * 70)

        for port_name in sorted(self._port_info.keys(), 
                                key=lambda x: int(x.replace('port', ''))):
            info = self._port_info[port_name]
            port_type = info['type'].upper()
            pos = info['position']
            n_faces = info['num_faces']

            type_marker = "◀▶" if info['type'] == 'external' else "◀┃▶"

            components = info.get('connected_components', [])
            comp_str = ", ".join(f"{k}.{p}" for k, p in components) if components else "?"

            print(f"  {port_name:8s} │ {self.main_axis}={pos:+.6f} │ "
                  f"{port_type:10s} │ {n_faces} face(s) {type_marker}")
            print(f"           │ Components: {comp_str}")

        print("=" * 70)

        # Visual schematic
        self._print_schematic()

    def _print_schematic(self) -> None:
        """Print a visual schematic of the assembly."""
        ports_sorted = sorted(
            self._port_info.items(), 
            key=lambda x: x[1]['position']
        )

        print(f"\nSchematic (along {self.main_axis}):\n")

        # Build port line
        port_parts = []
        for port_name, info in ports_sorted:
            if info['type'] == 'external':
                port_parts.append(f"║{port_name}║")
            else:
                port_parts.append(f"┃{port_name}┃")

        # Build component line
        comp_parts = []
        for key in self._component_order:
            entry = self._components[key]
            if entry.base_name != key:
                # Show base name for duplicates
                comp_parts.append(f"[{key}]")
            else:
                comp_parts.append(f"[{key}]")

        print("  Ports:      " + " ─── ".join(port_parts))
        print("  Components: " + " ─── ".join(comp_parts))

        # Show identical component groups
        identical = self.get_identical_components()
        if identical:
            print(f"\n  Identical components (compute once, reuse):")
            for base_name, keys in identical.items():
                print(f"    {base_name}: {', '.join(keys)}")
        print()

    # =========================================================================
    # Visualization
    # =========================================================================

    def inspect(
        self,
        show_labels: bool = True,
        color_by_component: bool = True,
        **kwargs
    ) -> None:
        """Visualize assembly configuration."""
        self.compute_layout()

        if not self._components:
            print("Assembly is empty")
            return

        shapes = []
        colors = self._generate_colors(len(self._components))

        for i, key in enumerate(self._component_order):
            entry = self._components[key]
            shape = entry.geometry.geo
            transformed = self._apply_transform(shape, entry.transform)

            if color_by_component:
                try:
                    for face in transformed.faces:
                        face.col = colors[i]
                except:
                    pass

            shapes.append(transformed)

        if len(shapes) == 1:
            display_geo = shapes[0]
        else:
            display_geo = Glue(shapes)

        Draw(display_geo, **kwargs)

        # Summary
        print(f"\nAssembly: {len(self._components)} components")
        print("-" * 60)
        for i, key in enumerate(self._component_order):
            entry = self._components[key]
            pos = entry.transform.translation
            geo_type = type(entry.geometry).__name__
            size = entry.size

            base_info = f" (={entry.base_name})" if entry.base_name != key else ""
            sub = " [sub-assembly]" if entry.is_assembly else ""
            cached = " [CACHED]" if entry.geometry.has_cached_solution() else ""

            print(f"  [{i}] {key}{base_info}: {geo_type}{sub}{cached}")
            print(f"      Position: ({pos[0]:.4f}, {pos[1]:.4f}, {pos[2]:.4f})")
            print(f"      Size:     ({size[0]:.4f}, {size[1]:.4f}, {size[2]:.4f})")

        # Show identical groups
        identical = self.get_identical_components()
        if identical:
            print(f"\nIdentical components:")
            for base_name, keys in identical.items():
                print(f"  {base_name}: {len(keys)} instances -> compute once")

    def show(
        self,
        what: Literal["geometry", "mesh", "geo", "inspect"] = "geometry",
        **kwargs
    ) -> None:
        """Display the assembly."""
        what = what.lower()

        if what == "inspect":
            self.inspect(**kwargs)
            return
        elif what in ("geometry", "geo"):
            if self.geo is None:
                raise ValueError("Assembly not built. Call build() first.")
            scene = Draw(self.geo, **kwargs)
        elif what == "mesh":
            if self.mesh is None:
                raise ValueError("Mesh not generated. Call generate_mesh() first.")
            scene = Draw(self.mesh, **kwargs)
        else:
            raise ValueError(f"Invalid option '{what}'")

        _display_webgui_fallback(scene)

    @staticmethod
    def _generate_colors(n: int) -> List[Tuple[float, float, float]]:
        """Generate distinct colors."""
        colors = []
        for i in range(n):
            hue = i / max(n, 1)
            h = hue * 6
            c = 0.9 * 0.7
            x = c * (1 - abs(h % 2 - 1))
            m = 0.9 - c

            if h < 1: r, g, b = c, x, 0
            elif h < 2: r, g, b = x, c, 0
            elif h < 3: r, g, b = 0, c, x
            elif h < 4: r, g, b = 0, x, c
            elif h < 5: r, g, b = x, 0, c
            else: r, g, b = c, 0, x

            colors.append((r + m, g + m, b + m))
        return colors

    # =========================================================================
    # Assembly Bounds and Info
    # =========================================================================

    def get_assembly_bounds(self) -> Tuple[Tuple[float, ...], Tuple[float, ...]]:
        """Get bounding box of entire assembly."""
        if not self._components:
            return ((0, 0, 0), (0, 0, 0))

        all_min = [float('inf')] * 3
        all_max = [float('-inf')] * 3

        for key in self._components:
            entry = self._components[key]
            pmin, pmax = entry.original_bounds or ((0, 0, 0), (1, 1, 1))
            t = entry.transform.translation

            for i in range(3):
                all_min[i] = min(all_min[i], pmin[i] + t[i])
                all_max[i] = max(all_max[i], pmax[i] + t[i])

        return tuple(all_min), tuple(all_max)

    def summary(self) -> Dict[str, Any]:
        """Get assembly summary information."""
        identical = self.get_identical_components()

        return {
            'n_components': len(self._components),
            'n_unique_geometries': len(self._base_name_groups),
            'n_identical_groups': len(identical),
            'n_connections': len(self._connections),
            'n_ports': len(self._port_info),
            'n_external_ports': len(self.get_external_ports()),
            'n_interface_ports': len(self.get_interface_ports()),
            'main_axis': self.main_axis,
            'is_built': self._is_built,
            'has_mesh': self.mesh is not None,
            'component_keys': list(self._component_order),
            'identical_groups': identical,
            'bounds': self.get_assembly_bounds() if self._components else None,
        }

    def print_info(self) -> None:
        """Print detailed assembly information."""
        info = self.summary()

        print("\n" + "=" * 70)
        print("ASSEMBLY INFORMATION")
        print("=" * 70)
        print(f"Total components:       {info['n_components']}")
        print(f"Unique geometries:      {info['n_unique_geometries']}")
        print(f"Identical groups:       {info['n_identical_groups']}")
        print(f"Main axis:              {info['main_axis']}")
        print(f"Built:                  {info['is_built']}")
        print(f"Has mesh:               {info['has_mesh']}")

        print(f"\nPorts:")
        print(f"  Total:                {info['n_ports']}")
        print(f"  External:             {info['n_external_ports']}")
        print(f"  Interfaces:           {info['n_interface_ports']}")

        if info['bounds']:
            pmin, pmax = info['bounds']
            print(f"\nBounding Box:")
            print(f"  Min: ({pmin[0]:.4f}, {pmin[1]:.4f}, {pmin[2]:.4f})")
            print(f"  Max: ({pmax[0]:.4f}, {pmax[1]:.4f}, {pmax[2]:.4f})")

        print(f"\nComponents:")
        for i, key in enumerate(self._component_order):
            entry = self._components[key]
            geo_type = type(entry.geometry).__name__
            pos = entry.transform.translation

            base_info = f" (base: {entry.base_name})" if entry.base_name != key else ""

            print(f"  [{i}] {key}{base_info}: {geo_type}")
            print(f"       Position: ({pos[0]:.4f}, {pos[1]:.4f}, {pos[2]:.4f})")

        # Identical component groups
        if info['identical_groups']:
            print(f"\nIdentical Component Groups (solver optimization):")
            for base_name, keys in info['identical_groups'].items():
                print(f"  '{base_name}': {len(keys)} instances")
                print(f"    Keys: {', '.join(keys)}")
                print(f"    -> Compute once, reuse {len(keys)-1} times")

        print("=" * 70)

        # Print port map if built
        if self._is_built and self._port_info:
            self.print_port_info()

    # =========================================================================
    # Tagging and Cache Support
    # =========================================================================

    def _get_geometry_params(self) -> Dict[str, Any]:
        """Get parameters for tag computation."""
        component_tags = {}
        for key, entry in self._components.items():
            component_tags[key] = {
                'tag_hash': entry.tag.full_hash,
                'base_name': entry.base_name,
                'transform': {
                    'translation': entry.transform.translation,
                    'rotation': entry.transform.rotation
                }
            }

        return {
            'class': 'Assembly',
            'main_axis': self.main_axis,
            'components': component_tags,
            'order': list(self._component_order),
            'identical_groups': self.get_identical_components()
        }

    def _get_mesh_params(self) -> Dict[str, Any]:
        """Get mesh parameters for tag computation."""
        return {
            'maxh': getattr(self, 'maxh', None),
            'component_count': len(self._components)
        }

    def get_solver_optimization_info(self) -> Dict[str, Any]:
        """
        Get information for solver optimization.

        Returns details about which components are identical and can
        share computed solutions.

        Returns
        -------
        dict
            {
                'total_components': int,
                'unique_computations': int,
                'reusable_results': int,
                'groups': {
                    base_name: {
                        'keys': [key1, key2, ...],
                        'compute_key': key1,  # Which one to compute
                        'reuse_keys': [key2, ...]  # Which ones reuse the result
                    }
                }
            }
        """
        identical = self.get_identical_components()

        # Single-instance components
        unique_keys = [
            keys[0] for name, keys in self._base_name_groups.items()
            if len(keys) == 1
        ]

        groups = {}
        total_reuse = 0

        for base_name, keys in identical.items():
            groups[base_name] = {
                'keys': keys,
                'compute_key': keys[0],
                'reuse_keys': keys[1:]
            }
            total_reuse += len(keys) - 1

        return {
            'total_components': len(self._components),
            'unique_computations': len(self._base_name_groups),
            'reusable_results': total_reuse,
            'efficiency': 1 - (len(self._base_name_groups) / len(self._components)) 
                         if self._components else 0,
            'unique_keys': unique_keys,
            'groups': groups
        }

    def save_geometry(self, project_path) -> None:
        """
        Save assembly geometry and all sub-component source files.

        Creates:
        - ``geometry/components/{key}_source.{ext}`` for imported CAD files
        - ``geometry/components/{key}.step`` for primitives
        - ``geometry/history.json`` with rewritten paths
        """
        project_path = Path(project_path)
        geo_dir = project_path / 'geometry'
        comp_dir = geo_dir / 'components'
        geo_dir.mkdir(parents=True, exist_ok=True)
        comp_dir.mkdir(parents=True, exist_ok=True)

        # Also export the whole assembly as a STEP file for external tools
        if self.geo is not None:
            try:
                self.geo.WriteStep(str(geo_dir / 'assembly.step'))
            except Exception:
                pass

        # Build history with rewritten sub-component paths
        history_for_save = []
        component_sources = {}  # key -> {source_link, source_hash, source_filename}

        for entry in self._history:
            e = dict(entry)

            if e.get('op') == 'add' or e.get('op') == 'connect':
                name = e.get('name', '')
                geo_type = e.get('geometry_type', '')
                sub_history = e.get('geometry_history', [])

                # Find the component key for this name
                # (might be suffixed like name_1, name_2)
                comp_key = None
                for k, comp in self._components.items():
                    if comp.base_name == name:
                        if k not in component_sources:
                            comp_key = k
                            break
                if comp_key is None:
                    comp_key = name

                comp_entry = self._components.get(comp_key)
                if comp_entry is not None:
                    geo_obj = comp_entry.geometry
                    from .importers import OCCImporter

                    if isinstance(geo_obj, OCCImporter) and geo_obj._source_link:
                        # Copy source CAD file
                        src = Path(geo_obj._source_link)
                        if src.exists():
                            ext = src.suffix
                            local_name = f'{comp_key}_source{ext}'
                            dest = comp_dir / local_name
                            if src.absolute() != dest.absolute():
                                shutil.copy2(str(src), str(dest))
                            source_hash = self._file_hash(dest)
                            component_sources[comp_key] = {
                                'source_link': str(geo_obj._source_link),
                                'source_hash': source_hash,
                                'source_filename': f'components/{local_name}',
                            }
                            # Rewrite filepath in sub-history
                            rewritten_sub = []
                            for sh in sub_history:
                                sh_copy = dict(sh)
                                if sh_copy.get('op') in ('import_occ', 'import_step'):
                                    sh_copy['filepath'] = f'geometry/components/{local_name}'
                                rewritten_sub.append(sh_copy)
                            e['geometry_history'] = rewritten_sub

                    elif not isinstance(geo_obj, Assembly):
                        # Primitive — export as STEP
                        if geo_obj.geo is not None:
                            try:
                                step_name = f'{comp_key}.step'
                                geo_obj.geo.WriteStep(str(comp_dir / step_name))
                                component_sources[comp_key] = {
                                    'source_link': None,
                                    'source_hash': None,
                                    'source_filename': f'components/{step_name}',
                                }
                            except Exception:
                                pass

            history_for_save.append(e)

        meta = {
            'type': 'Assembly',
            'module': 'cavsim3d.geometry.assembly',
            'source_link': None,
            'source_filename': None,
            'source_hash': None,
            'component_sources': component_sources,
            'history': history_for_save,
        }

        with open(geo_dir / 'history.json', 'w') as f:
            json.dump(meta, f, indent=2, default=str)

    @classmethod
    def _rebuild_from_history(
        cls,
        history: List[dict],
        project_path,
        source_file=None,
    ) -> 'Assembly':
        """Reconstruct assembly by replaying operation history."""
        project_path = Path(project_path)
        obj = None

        for entry in history:
            op = entry['op']
            params = {k: v for k, v in entry.items() if k not in ['op', 'timestamp']}

            if op == '__init__':
                obj = cls(main_axis=params.get('main_axis', 'Z'))

            elif op == 'add':
                # Reconstruct sub-geometry
                sub_history = list(params.pop('geometry_history'))
                sub_type = params.pop('geometry_type')

                sub_cls = BaseGeometry._get_subclass(sub_type)
                if sub_cls is None:
                    raise ValueError(f"Unknown geometry type '{sub_type}'")

                # Resolve source file from rewritten paths
                sub_source = None
                for sh in sub_history:
                    if sh.get('op') in ('import_occ', 'import_step'):
                        fp = sh.get('filepath', '')
                        candidate = project_path / fp
                        if candidate.exists():
                            sub_source = candidate
                        break

                sub_geo = sub_cls._rebuild_from_history(
                    sub_history, project_path, source_file=sub_source
                )

                name = params.pop('name')
                obj.add(name, sub_geo, **params)

            elif op == 'connect':
                sub_history = list(params.pop('geometry_history'))
                sub_type = params.pop('geometry_type')

                sub_cls = BaseGeometry._get_subclass(sub_type)
                if sub_cls is None:
                    raise ValueError(f"Unknown geometry type '{sub_type}'")

                sub_source = None
                for sh in sub_history:
                    if sh.get('op') in ('import_occ', 'import_step'):
                        fp = sh.get('filepath', '')
                        candidate = project_path / fp
                        if candidate.exists():
                            sub_source = candidate
                        break

                sub_geo = sub_cls._rebuild_from_history(
                    sub_history, project_path, source_file=sub_source
                )

                name = params.pop('name')
                obj.connect(name, sub_geo, **params)

            elif op == 'remove':
                obj.remove(**params)

            elif op == 'rotate':
                obj.rotate(**params)

            elif op == 'translate':
                obj.translate(**params)

            elif op == 'build':
                obj.build()

            elif op == 'generate_mesh':
                obj.generate_mesh(**params)

        return obj

Functions

__init__(main_axis='Z')

Source code in cavsim3d/geometry/assembly.py

def __init__(self, main_axis: Literal['X', 'Y', 'Z'] = 'Z'):
    super().__init__()
    self.main_axis = main_axis.upper()

    self._components: Dict[str, ComponentEntry] = {}
    self._component_order: List[str] = []
    self._connections: List[Connection] = []

    # Track base names for duplicate detection
    self._base_name_counts: Dict[str, int] = {}
    self._base_name_groups: Dict[str, List[str]] = {}  # base_name -> [key1, key2, ...]

    self._layout_computed = False
    self._is_built = False

    # Port and solid information (populated by build())
    self._port_info: Dict[str, Dict] = {}
    self._solid_info: Dict[str, Dict] = {}

    self._record('__init__', main_axis=main_axis)

add(name, geometry, position=None, rotation=None, after=None, before=None, align_port=None, **metadata)

Add a component (geometry or sub-assembly) to the assembly.

Duplicate names are allowed and automatically suffixed with _1, _2, etc. Components with the same base name are tracked for solver optimization (compute once, reuse for identical geometries).

Parameters

name : str Component name (duplicates allowed, will be auto-suffixed) geometry : BaseGeometry or Assembly Component to add (can be another Assembly!) position : tuple, optional Explicit position (x, y, z) rotation : tuple, optional Rotation angles in degrees (rx, ry, rz) after : str, optional Place after this component along main axis (can use base name) before : str, optional Place before this component along main axis (can use base name) align_port : str, optional Port to use for alignment **metadata Additional metadata

Returns

self : Assembly

Examples

assembly.add("midcell", midcell) # -> "midcell_1" assembly.add("midcell", midcell) # -> "midcell_2"
assembly.add("midcell", midcell) # -> "midcell_3" assembly.get_identical_components()

Source code in cavsim3d/geometry/assembly.py

def add(
    self,
    name: str,
    geometry: Union[BaseGeometry, 'Assembly'],
    position: Optional[Tuple[float, float, float]] = None,
    rotation: Optional[Tuple[float, float, float]] = None,
    after: Optional[str] = None,
    before: Optional[str] = None,
    align_port: Optional[str] = None,
    **metadata
) -> 'Assembly':
    """
    Add a component (geometry or sub-assembly) to the assembly.

    Duplicate names are allowed and automatically suffixed with _1, _2, etc.
    Components with the same base name are tracked for solver optimization
    (compute once, reuse for identical geometries).

    Parameters
    ----------
    name : str
        Component name (duplicates allowed, will be auto-suffixed)
    geometry : BaseGeometry or Assembly
        Component to add (can be another Assembly!)
    position : tuple, optional
        Explicit position (x, y, z)
    rotation : tuple, optional
        Rotation angles in degrees (rx, ry, rz)
    after : str, optional
        Place after this component along main axis (can use base name)
    before : str, optional
        Place before this component along main axis (can use base name)
    align_port : str, optional
        Port to use for alignment
    **metadata
        Additional metadata

    Returns
    -------
    self : Assembly

    Examples
    --------
    >>> assembly.add("midcell", midcell)           # -> "midcell_1"
    >>> assembly.add("midcell", midcell)           # -> "midcell_2"  
    >>> assembly.add("midcell", midcell)           # -> "midcell_3"
    >>> assembly.get_identical_components()
    {'midcell': ['midcell_1', 'midcell_2', 'midcell_3']}
    """
    base_name = name

    # Generate unique key if name already exists
    if name in self._components or name in self._base_name_counts:
        # This is a duplicate name
        if base_name not in self._base_name_counts:
            # First duplicate - rename the original
            self._rename_first_instance(base_name)

        # Increment count and generate new key
        self._base_name_counts[base_name] += 1
        key = f"{base_name}_{self._base_name_counts[base_name]}"
    else:
        # First instance of this name
        key = name
        self._base_name_counts[base_name] = 0  # Will become 1 if duplicated

    # Track in base name groups
    if base_name not in self._base_name_groups:
        self._base_name_groups[base_name] = []
    self._base_name_groups[base_name].append(key)

    # Handle sub-assemblies
    if isinstance(geometry, Assembly):
        if not geometry._is_built:
            geometry.build()
    elif geometry.geo is None:
        # Auto-build if possible instead of raising error
        try:
            geometry.build()
        except Exception as e:
            raise ValueError(f"Geometry '{name}' not built and auto-build failed: {e}. Call build() first.")

    transform = Transform3D(
        translation=position or (0.0, 0.0, 0.0),
        rotation=rotation or (0.0, 0.0, 0.0)
    )

    entry = ComponentEntry(
        geometry=geometry,
        key=key,
        base_name=base_name,
        transform=transform,
        aligned_port=align_port,
        metadata=metadata
    )

    self._components[key] = entry

    # Handle ordering - resolve component and optional gap
    gap = 0.0
    if isinstance(after, tuple):
        ref_after, gap = after
    else:
        ref_after = after

    if isinstance(before, tuple):
        ref_before, gap = before
    else:
        ref_before = before

    resolved_after = self._resolve_component_ref(ref_after) if ref_after else None
    resolved_before = self._resolve_component_ref(ref_before) if ref_before else None

    if resolved_after is not None:
        idx = self._component_order.index(resolved_after) + 1
        self._component_order.insert(idx, key)

        self._connections.append(Connection(
            from_key=resolved_after,
            to_key=key,
            from_port='port2',
            to_port=align_port or 'port1',
            gap=gap
        ))

    elif resolved_before is not None:
        idx = self._component_order.index(resolved_before)
        self._component_order.insert(idx, key)

        self._connections.append(Connection(
            from_key=key,
            to_key=resolved_before,
            from_port=align_port or 'port2',
            to_port='port1',
            gap=-gap  # Before means moving current component backwards relative to target
        ))
    else:
        self._component_order.append(key)

    self._layout_computed = False
    self._is_built = False
    self.invalidate_tag()
    self._record(
        'add', 
        name=name, 
        geometry_type=type(geometry).__name__,
        geometry_history=geometry.get_history(),
        position=position,
        rotation=rotation,
        after=after,
        before=before,
        align_port=align_port,
        **metadata
    )

    return self

remove(ref)

Remove a component by key or base name (removes all instances if base name).

Source code in cavsim3d/geometry/assembly.py

def remove(self, ref: str) -> 'Assembly':
    """Remove a component by key or base name (removes all instances if base name)."""
    keys_to_remove = []

    if ref in self._components:
        keys_to_remove = [ref]
    elif ref in self._base_name_groups:
        keys_to_remove = list(self._base_name_groups[ref])
    else:
        raise KeyError(f"Component '{ref}' not found")

    for key in keys_to_remove:
        entry = self._components[key]
        base_name = entry.base_name

        del self._components[key]
        self._component_order.remove(key)
        self._connections = [
            c for c in self._connections 
            if c.from_key != key and c.to_key != key
        ]

        # Update base name tracking
        if base_name in self._base_name_groups:
            self._base_name_groups[base_name].remove(key)
            if not self._base_name_groups[base_name]:
                del self._base_name_groups[base_name]
                del self._base_name_counts[base_name]

    self._layout_computed = False
    self._is_built = False
    self.invalidate_tag()

    self._record('remove', ref=ref)

    return self

connect(name, geometry, to_component, from_port='port1', to_port='port2', rotation=None, gap=0.0, **metadata)

Add component connected via ports to existing component.

Source code in cavsim3d/geometry/assembly.py

def connect(
    self,
    name: str,
    geometry: Union[BaseGeometry, 'Assembly'],
    to_component: str,
    from_port: str = "port1",
    to_port: str = "port2",
    rotation: Optional[Tuple[float, float, float]] = None,
    gap: float = 0.0,
    **metadata
) -> 'Assembly':
    """Add component connected via ports to existing component."""
    resolved_to = self._resolve_component_ref(to_component)

    self.add(name, geometry, rotation=rotation, **metadata)

    # Get the key that was just added (might be suffixed)
    new_key = self._component_order[-1]

    self._record(
        'connect',
        name=name,
        geometry_type=type(geometry).__name__,
        geometry_history=geometry.get_history(),
        to_component=to_component,
        from_port=from_port,
        to_port=to_port,
        rotation=rotation,
        gap=gap,
        **metadata
    )

    return self

build()

Build multi-solid geometry from all components.

Creates a compound where each component is a separate solid with: - Solid named after component key (e.g., 'midcell_1', 'midcell_2') - Sequential port naming (port1, port2, ...) along main axis - Interface ports (where solids meet) shared between adjacent components - External ports at assembly ends - All other faces named 'wall'

Source code in cavsim3d/geometry/assembly.py

def build(self) -> None:
    """
    Build multi-solid geometry from all components.

    Creates a compound where each component is a separate solid with:
    - Solid named after component key (e.g., 'midcell_1', 'midcell_2')
    - Sequential port naming (port1, port2, ...) along main axis
    - Interface ports (where solids meet) shared between adjacent components
    - External ports at assembly ends
    - All other faces named 'wall'
    """
    if not self._components:
        raise ValueError("No components in assembly")

    self.compute_layout()

    # NEW: Snap interfaces to compensate for CAD inaccuracies
    self._snap_interfaces()

    shapes = []
    # Build sub-assemblies first
    for entry in self._components.values():
        if entry.is_assembly and not entry.geometry._is_built:
            entry.geometry.build()

    # Collect transformed shapes
    shapes = []
    shape_keys = []  # Track which key each shape belongs to

    for key in self._component_order:
        entry = self._components[key]
        shape = entry.geometry.geo
        transformed = self._apply_transform(shape, entry.transform)
        shapes.append(transformed)
        shape_keys.append(key)

    # To prevent merging identical/co-planar solids, we assign unique 
    # temporary materials and use OCCGeometry(shapes).shape which is
    # often more robust than Glue() for preserving solid identities.
    if len(shapes) > 1:
        for i, shape in enumerate(shapes):
            shape.mat(f"__temp_solid_{i}")

    # Create multi-solid compound
    if len(shapes) == 1:
        self.geo = shapes[0]
    else:
        # OCCGeometry(list) performs gluing but preserves solid domains better
        occ_geo = OCCGeometry(shapes)
        self.geo = occ_geo.shape

    # Name solids and setup boundaries
    self._name_solids(shape_keys)
    self._setup_boundaries()

    self._is_built = True
    self._record('build')

generate_mesh(maxh=None, curve_order=3)

Generate mesh with support for per-component refinement.

Source code in cavsim3d/geometry/assembly.py

def generate_mesh(self, maxh=None, curve_order: int = 3) -> Mesh:
    """
    Generate mesh with support for per-component refinement.
    """
    if self.geo is None:
        self.build()

    # Apply per-component maxh if available
    solids = list(self.geo.solids)
    if len(solids) == len(self._component_order):
        for i, key in enumerate(self._component_order):
            entry = self._components[key]
            # Try to get maxh from component geometry
            comp_maxh = getattr(entry.geometry, 'maxh', None)
            if comp_maxh:
                solids[i].maxh = comp_maxh

    # Call base generate_mesh which will create OCCGeometry(self.geo)
    # and respect the per-solid maxh settings.
    return super().generate_mesh(maxh=maxh, curve_order=curve_order)

compute_layout()

Compute positions based on connections.

Source code in cavsim3d/geometry/assembly.py

def compute_layout(self) -> 'Assembly':
    """Compute positions based on connections."""
    if self._layout_computed:
        return self

    if not self._components:
        self._layout_computed = True
        return self

    positioned = set()
    first_key = self._component_order[0]
    positioned.add(first_key)

    for conn in self._connections:
        if conn.to_key in positioned and conn.from_key in positioned:
            continue
        if conn.from_key in positioned:
            self._position_connected(conn)
            positioned.add(conn.to_key)
        elif conn.to_key in positioned:
            # Reverse connection positioning not yet implemented for all cases
            pass

    # Safety for unpositioned components
    for key in self._component_order:
        if key not in positioned:
            if positioned:
                last = list(positioned)[-1]
                self._position_after(key, last)
            positioned.add(key)

    self._layout_computed = True
    return self

rotate(ref, angles, center=None)

Apply rotation to a component.

Source code in cavsim3d/geometry/assembly.py

def rotate(
    self, 
    ref: str, 
    angles: Tuple[float, float, float],
    center: Optional[Tuple[float, float, float]] = None
) -> 'Assembly':
    """Apply rotation to a component."""
    key = self._resolve_component_ref(ref)
    entry = self._components[key]
    current = entry.transform.rotation
    entry.transform = Transform3D(
        translation=entry.transform.translation,
        rotation=tuple(c + a for c, a in zip(current, angles)),
        rotation_center=center
    )
    self._layout_computed = False
    self._is_built = False
    return self

translate(ref, offset)

Apply translation to a component.

Source code in cavsim3d/geometry/assembly.py

def translate(self, ref: str, offset: Tuple[float, float, float]) -> 'Assembly':
    """Apply translation to a component."""
    key = self._resolve_component_ref(ref)
    entry = self._components[key]
    current = entry.transform.translation
    entry.transform = Transform3D(
        translation=tuple(c + o for c, o in zip(current, offset)),
        rotation=entry.transform.rotation,
        rotation_center=entry.transform.rotation_center
    )
    self._layout_computed = False
    self._is_built = False
    return self

show(what='geometry', **kwargs)

Display the assembly.

Source code in cavsim3d/geometry/assembly.py

def show(
    self,
    what: Literal["geometry", "mesh", "geo", "inspect"] = "geometry",
    **kwargs
) -> None:
    """Display the assembly."""
    what = what.lower()

    if what == "inspect":
        self.inspect(**kwargs)
        return
    elif what in ("geometry", "geo"):
        if self.geo is None:
            raise ValueError("Assembly not built. Call build() first.")
        scene = Draw(self.geo, **kwargs)
    elif what == "mesh":
        if self.mesh is None:
            raise ValueError("Mesh not generated. Call generate_mesh() first.")
        scene = Draw(self.mesh, **kwargs)
    else:
        raise ValueError(f"Invalid option '{what}'")

    _display_webgui_fallback(scene)

inspect(show_labels=True, color_by_component=True, **kwargs)

Visualize assembly configuration.

Source code in cavsim3d/geometry/assembly.py

def inspect(
    self,
    show_labels: bool = True,
    color_by_component: bool = True,
    **kwargs
) -> None:
    """Visualize assembly configuration."""
    self.compute_layout()

    if not self._components:
        print("Assembly is empty")
        return

    shapes = []
    colors = self._generate_colors(len(self._components))

    for i, key in enumerate(self._component_order):
        entry = self._components[key]
        shape = entry.geometry.geo
        transformed = self._apply_transform(shape, entry.transform)

        if color_by_component:
            try:
                for face in transformed.faces:
                    face.col = colors[i]
            except:
                pass

        shapes.append(transformed)

    if len(shapes) == 1:
        display_geo = shapes[0]
    else:
        display_geo = Glue(shapes)

    Draw(display_geo, **kwargs)

    # Summary
    print(f"\nAssembly: {len(self._components)} components")
    print("-" * 60)
    for i, key in enumerate(self._component_order):
        entry = self._components[key]
        pos = entry.transform.translation
        geo_type = type(entry.geometry).__name__
        size = entry.size

        base_info = f" (={entry.base_name})" if entry.base_name != key else ""
        sub = " [sub-assembly]" if entry.is_assembly else ""
        cached = " [CACHED]" if entry.geometry.has_cached_solution() else ""

        print(f"  [{i}] {key}{base_info}: {geo_type}{sub}{cached}")
        print(f"      Position: ({pos[0]:.4f}, {pos[1]:.4f}, {pos[2]:.4f})")
        print(f"      Size:     ({size[0]:.4f}, {size[1]:.4f}, {size[2]:.4f})")

    # Show identical groups
    identical = self.get_identical_components()
    if identical:
        print(f"\nIdentical components:")
        for base_name, keys in identical.items():
            print(f"  {base_name}: {len(keys)} instances -> compute once")

get_port_info()

Get information about all ports in the assembly.

Source code in cavsim3d/geometry/assembly.py

def get_port_info(self) -> Dict[str, Dict]:
    """Get information about all ports in the assembly."""
    return dict(self._port_info) if self._port_info else {}

get_solid_info()

Get information about all solids in the assembly.

Source code in cavsim3d/geometry/assembly.py

def get_solid_info(self) -> Dict[str, Dict]:
    """Get information about all solids in the assembly."""
    return dict(self._solid_info) if self._solid_info else {}

get_external_ports()

Get list of external port names.

Source code in cavsim3d/geometry/assembly.py

def get_external_ports(self) -> List[str]:
    """Get list of external port names."""
    return [name for name, info in self._port_info.items() 
            if info['type'] == 'external']

get_interface_ports()

Get list of interface port names.

Source code in cavsim3d/geometry/assembly.py

def get_interface_ports(self) -> List[str]:
    """Get list of interface port names."""
    return [name for name, info in self._port_info.items() 
            if info['type'] == 'interface']

get_assembly_bounds()

Get bounding box of entire assembly.

Source code in cavsim3d/geometry/assembly.py

def get_assembly_bounds(self) -> Tuple[Tuple[float, ...], Tuple[float, ...]]:
    """Get bounding box of entire assembly."""
    if not self._components:
        return ((0, 0, 0), (0, 0, 0))

    all_min = [float('inf')] * 3
    all_max = [float('-inf')] * 3

    for key in self._components:
        entry = self._components[key]
        pmin, pmax = entry.original_bounds or ((0, 0, 0), (1, 1, 1))
        t = entry.transform.translation

        for i in range(3):
            all_min[i] = min(all_min[i], pmin[i] + t[i])
            all_max[i] = max(all_max[i], pmax[i] + t[i])

    return tuple(all_min), tuple(all_max)

get_identical_components()

Get groups of components that share the same base name (geometry).

These components are candidates for solver optimization - compute the solution once and reuse for all identical instances.

Returns

dict {base_name: [key1, key2, ...]} for groups with more than one member

Source code in cavsim3d/geometry/assembly.py

def get_identical_components(self) -> Dict[str, List[str]]:
    """
    Get groups of components that share the same base name (geometry).

    These components are candidates for solver optimization - compute
    the solution once and reuse for all identical instances.

    Returns
    -------
    dict
        {base_name: [key1, key2, ...]} for groups with more than one member
    """
    return {
        name: keys for name, keys in self._base_name_groups.items()
        if len(keys) > 1
    }

get_components_by_base_name(base_name)

Get all component keys with a given base name.

Source code in cavsim3d/geometry/assembly.py

def get_components_by_base_name(self, base_name: str) -> List[str]:
    """Get all component keys with a given base name."""
    return self._base_name_groups.get(base_name, [])

get_solver_optimization_info()

Get information for solver optimization.

Returns details about which components are identical and can share computed solutions.

Returns

dict { 'total_components': int, 'unique_computations': int, 'reusable_results': int, 'groups': { base_name: { 'keys': [key1, key2, ...], 'compute_key': key1, # Which one to compute 'reuse_keys': [key2, ...] # Which ones reuse the result } } }

Source code in cavsim3d/geometry/assembly.py

def get_solver_optimization_info(self) -> Dict[str, Any]:
    """
    Get information for solver optimization.

    Returns details about which components are identical and can
    share computed solutions.

    Returns
    -------
    dict
        {
            'total_components': int,
            'unique_computations': int,
            'reusable_results': int,
            'groups': {
                base_name: {
                    'keys': [key1, key2, ...],
                    'compute_key': key1,  # Which one to compute
                    'reuse_keys': [key2, ...]  # Which ones reuse the result
                }
            }
        }
    """
    identical = self.get_identical_components()

    # Single-instance components
    unique_keys = [
        keys[0] for name, keys in self._base_name_groups.items()
        if len(keys) == 1
    ]

    groups = {}
    total_reuse = 0

    for base_name, keys in identical.items():
        groups[base_name] = {
            'keys': keys,
            'compute_key': keys[0],
            'reuse_keys': keys[1:]
        }
        total_reuse += len(keys) - 1

    return {
        'total_components': len(self._components),
        'unique_computations': len(self._base_name_groups),
        'reusable_results': total_reuse,
        'efficiency': 1 - (len(self._base_name_groups) / len(self._components)) 
                     if self._components else 0,
        'unique_keys': unique_keys,
        'groups': groups
    }

print_port_info()

Print detailed port information.

Source code in cavsim3d/geometry/assembly.py

def print_port_info(self) -> None:
    """Print detailed port information."""
    if not self._port_info:
        print("No port information available. Call build() first.")
        return

    print("\n" + "=" * 70)
    print("ASSEMBLY PORT MAP")
    print("=" * 70)
    print(f"Main axis: {self.main_axis}")
    print(f"Total ports: {len(self._port_info)}")
    print("-" * 70)

    for port_name in sorted(self._port_info.keys(), 
                            key=lambda x: int(x.replace('port', ''))):
        info = self._port_info[port_name]
        port_type = info['type'].upper()
        pos = info['position']
        n_faces = info['num_faces']

        type_marker = "◀▶" if info['type'] == 'external' else "◀┃▶"

        components = info.get('connected_components', [])
        comp_str = ", ".join(f"{k}.{p}" for k, p in components) if components else "?"

        print(f"  {port_name:8s} │ {self.main_axis}={pos:+.6f} │ "
              f"{port_type:10s} │ {n_faces} face(s) {type_marker}")
        print(f"           │ Components: {comp_str}")

    print("=" * 70)

    # Visual schematic
    self._print_schematic()

print_info()

Print detailed assembly information.

Source code in cavsim3d/geometry/assembly.py

def print_info(self) -> None:
    """Print detailed assembly information."""
    info = self.summary()

    print("\n" + "=" * 70)
    print("ASSEMBLY INFORMATION")
    print("=" * 70)
    print(f"Total components:       {info['n_components']}")
    print(f"Unique geometries:      {info['n_unique_geometries']}")
    print(f"Identical groups:       {info['n_identical_groups']}")
    print(f"Main axis:              {info['main_axis']}")
    print(f"Built:                  {info['is_built']}")
    print(f"Has mesh:               {info['has_mesh']}")

    print(f"\nPorts:")
    print(f"  Total:                {info['n_ports']}")
    print(f"  External:             {info['n_external_ports']}")
    print(f"  Interfaces:           {info['n_interface_ports']}")

    if info['bounds']:
        pmin, pmax = info['bounds']
        print(f"\nBounding Box:")
        print(f"  Min: ({pmin[0]:.4f}, {pmin[1]:.4f}, {pmin[2]:.4f})")
        print(f"  Max: ({pmax[0]:.4f}, {pmax[1]:.4f}, {pmax[2]:.4f})")

    print(f"\nComponents:")
    for i, key in enumerate(self._component_order):
        entry = self._components[key]
        geo_type = type(entry.geometry).__name__
        pos = entry.transform.translation

        base_info = f" (base: {entry.base_name})" if entry.base_name != key else ""

        print(f"  [{i}] {key}{base_info}: {geo_type}")
        print(f"       Position: ({pos[0]:.4f}, {pos[1]:.4f}, {pos[2]:.4f})")

    # Identical component groups
    if info['identical_groups']:
        print(f"\nIdentical Component Groups (solver optimization):")
        for base_name, keys in info['identical_groups'].items():
            print(f"  '{base_name}': {len(keys)} instances")
            print(f"    Keys: {', '.join(keys)}")
            print(f"    -> Compute once, reuse {len(keys)-1} times")

    print("=" * 70)

    # Print port map if built
    if self._is_built and self._port_info:
        self.print_port_info()

summary()

Get assembly summary information.

Source code in cavsim3d/geometry/assembly.py

def summary(self) -> Dict[str, Any]:
    """Get assembly summary information."""
    identical = self.get_identical_components()

    return {
        'n_components': len(self._components),
        'n_unique_geometries': len(self._base_name_groups),
        'n_identical_groups': len(identical),
        'n_connections': len(self._connections),
        'n_ports': len(self._port_info),
        'n_external_ports': len(self.get_external_ports()),
        'n_interface_ports': len(self.get_interface_ports()),
        'main_axis': self.main_axis,
        'is_built': self._is_built,
        'has_mesh': self.mesh is not None,
        'component_keys': list(self._component_order),
        'identical_groups': identical,
        'bounds': self.get_assembly_bounds() if self._components else None,
    }

save_geometry(project_path)

Save assembly geometry and all sub-component source files.

Creates: - geometry/components/{key}_source.{ext} for imported CAD files - geometry/components/{key}.step for primitives - geometry/history.json with rewritten paths

Source code in cavsim3d/geometry/assembly.py

def save_geometry(self, project_path) -> None:
    """
    Save assembly geometry and all sub-component source files.

    Creates:
    - ``geometry/components/{key}_source.{ext}`` for imported CAD files
    - ``geometry/components/{key}.step`` for primitives
    - ``geometry/history.json`` with rewritten paths
    """
    project_path = Path(project_path)
    geo_dir = project_path / 'geometry'
    comp_dir = geo_dir / 'components'
    geo_dir.mkdir(parents=True, exist_ok=True)
    comp_dir.mkdir(parents=True, exist_ok=True)

    # Also export the whole assembly as a STEP file for external tools
    if self.geo is not None:
        try:
            self.geo.WriteStep(str(geo_dir / 'assembly.step'))
        except Exception:
            pass

    # Build history with rewritten sub-component paths
    history_for_save = []
    component_sources = {}  # key -> {source_link, source_hash, source_filename}

    for entry in self._history:
        e = dict(entry)

        if e.get('op') == 'add' or e.get('op') == 'connect':
            name = e.get('name', '')
            geo_type = e.get('geometry_type', '')
            sub_history = e.get('geometry_history', [])

            # Find the component key for this name
            # (might be suffixed like name_1, name_2)
            comp_key = None
            for k, comp in self._components.items():
                if comp.base_name == name:
                    if k not in component_sources:
                        comp_key = k
                        break
            if comp_key is None:
                comp_key = name

            comp_entry = self._components.get(comp_key)
            if comp_entry is not None:
                geo_obj = comp_entry.geometry
                from .importers import OCCImporter

                if isinstance(geo_obj, OCCImporter) and geo_obj._source_link:
                    # Copy source CAD file
                    src = Path(geo_obj._source_link)
                    if src.exists():
                        ext = src.suffix
                        local_name = f'{comp_key}_source{ext}'
                        dest = comp_dir / local_name
                        if src.absolute() != dest.absolute():
                            shutil.copy2(str(src), str(dest))
                        source_hash = self._file_hash(dest)
                        component_sources[comp_key] = {
                            'source_link': str(geo_obj._source_link),
                            'source_hash': source_hash,
                            'source_filename': f'components/{local_name}',
                        }
                        # Rewrite filepath in sub-history
                        rewritten_sub = []
                        for sh in sub_history:
                            sh_copy = dict(sh)
                            if sh_copy.get('op') in ('import_occ', 'import_step'):
                                sh_copy['filepath'] = f'geometry/components/{local_name}'
                            rewritten_sub.append(sh_copy)
                        e['geometry_history'] = rewritten_sub

                elif not isinstance(geo_obj, Assembly):
                    # Primitive — export as STEP
                    if geo_obj.geo is not None:
                        try:
                            step_name = f'{comp_key}.step'
                            geo_obj.geo.WriteStep(str(comp_dir / step_name))
                            component_sources[comp_key] = {
                                'source_link': None,
                                'source_hash': None,
                                'source_filename': f'components/{step_name}',
                            }
                        except Exception:
                            pass

        history_for_save.append(e)

    meta = {
        'type': 'Assembly',
        'module': 'cavsim3d.geometry.assembly',
        'source_link': None,
        'source_filename': None,
        'source_hash': None,
        'component_sources': component_sources,
        'history': history_for_save,
    }

    with open(geo_dir / 'history.json', 'w') as f:
        json.dump(meta, f, indent=2, default=str)

Properties

Bases: BaseGeometry

Assembly of geometry components with multi-solid output.

Creates a compound geometry where each component remains a separate solid, similar to split geometries. Supports automatic positioning, unified port naming, and tracks identical components for solver optimization.

Parameters

main_axis : str Primary axis for concatenation ('X', 'Y', or 'Z')

Examples

TESLA 9-cell cavity with identical midcells:

assembly = Assembly(main_axis='Z') assembly.add("endcell_l", endcell_l) for i in range(7): ... assembly.add("midcell", midcell, after="endcell_l" if i == 0 else "midcell") assembly.add("endcell_r", endcell_r, after="midcell") assembly.build()

Check what was created

assembly.print_info()

Shows: midcell_1, midcell_2, ..., midcell_7 (all identical)

Resulting structure: - Solids: endcell_l, midcell_1, midcell_2, ..., midcell_7, endcell_r - Ports: port1 (external), port2 (interface), ..., port10 (external)

Source code in cavsim3d/geometry/assembly.py

class Assembly(BaseGeometry):
    """
    Assembly of geometry components with multi-solid output.

    Creates a compound geometry where each component remains a separate solid,
    similar to split geometries. Supports automatic positioning, unified port
    naming, and tracks identical components for solver optimization.

    Parameters
    ----------
    main_axis : str
        Primary axis for concatenation ('X', 'Y', or 'Z')

    Examples
    --------
    **TESLA 9-cell cavity with identical midcells:**

    >>> assembly = Assembly(main_axis='Z')
    >>> assembly.add("endcell_l", endcell_l)
    >>> for i in range(7):
    ...     assembly.add("midcell", midcell, after="endcell_l" if i == 0 else "midcell")
    >>> assembly.add("endcell_r", endcell_r, after="midcell")
    >>> assembly.build()
    >>> 
    >>> # Check what was created
    >>> assembly.print_info()
    >>> # Shows: midcell_1, midcell_2, ..., midcell_7 (all identical)

    **Resulting structure:**
    - Solids: endcell_l, midcell_1, midcell_2, ..., midcell_7, endcell_r
    - Ports: port1 (external), port2 (interface), ..., port10 (external)
    """

    _AXIS_VEC = {'X': X, 'Y': Y, 'Z': Z}
    _AXIS_IDX = {'X': 0, 'Y': 1, 'Z': 2}

    def __init__(self, main_axis: Literal['X', 'Y', 'Z'] = 'Z'):
        super().__init__()
        self.main_axis = main_axis.upper()

        self._components: Dict[str, ComponentEntry] = {}
        self._component_order: List[str] = []
        self._connections: List[Connection] = []

        # Track base names for duplicate detection
        self._base_name_counts: Dict[str, int] = {}
        self._base_name_groups: Dict[str, List[str]] = {}  # base_name -> [key1, key2, ...]

        self._layout_computed = False
        self._is_built = False

        # Port and solid information (populated by build())
        self._port_info: Dict[str, Dict] = {}
        self._solid_info: Dict[str, Dict] = {}

        self._record('__init__', main_axis=main_axis)

    # =========================================================================
    # Component Management (with duplicate name support)
    # =========================================================================

    def add(
        self,
        name: str,
        geometry: Union[BaseGeometry, 'Assembly'],
        position: Optional[Tuple[float, float, float]] = None,
        rotation: Optional[Tuple[float, float, float]] = None,
        after: Optional[str] = None,
        before: Optional[str] = None,
        align_port: Optional[str] = None,
        **metadata
    ) -> 'Assembly':
        """
        Add a component (geometry or sub-assembly) to the assembly.

        Duplicate names are allowed and automatically suffixed with _1, _2, etc.
        Components with the same base name are tracked for solver optimization
        (compute once, reuse for identical geometries).

        Parameters
        ----------
        name : str
            Component name (duplicates allowed, will be auto-suffixed)
        geometry : BaseGeometry or Assembly
            Component to add (can be another Assembly!)
        position : tuple, optional
            Explicit position (x, y, z)
        rotation : tuple, optional
            Rotation angles in degrees (rx, ry, rz)
        after : str, optional
            Place after this component along main axis (can use base name)
        before : str, optional
            Place before this component along main axis (can use base name)
        align_port : str, optional
            Port to use for alignment
        **metadata
            Additional metadata

        Returns
        -------
        self : Assembly

        Examples
        --------
        >>> assembly.add("midcell", midcell)           # -> "midcell_1"
        >>> assembly.add("midcell", midcell)           # -> "midcell_2"  
        >>> assembly.add("midcell", midcell)           # -> "midcell_3"
        >>> assembly.get_identical_components()
        {'midcell': ['midcell_1', 'midcell_2', 'midcell_3']}
        """
        base_name = name

        # Generate unique key if name already exists
        if name in self._components or name in self._base_name_counts:
            # This is a duplicate name
            if base_name not in self._base_name_counts:
                # First duplicate - rename the original
                self._rename_first_instance(base_name)

            # Increment count and generate new key
            self._base_name_counts[base_name] += 1
            key = f"{base_name}_{self._base_name_counts[base_name]}"
        else:
            # First instance of this name
            key = name
            self._base_name_counts[base_name] = 0  # Will become 1 if duplicated

        # Track in base name groups
        if base_name not in self._base_name_groups:
            self._base_name_groups[base_name] = []
        self._base_name_groups[base_name].append(key)

        # Handle sub-assemblies
        if isinstance(geometry, Assembly):
            if not geometry._is_built:
                geometry.build()
        elif geometry.geo is None:
            # Auto-build if possible instead of raising error
            try:
                geometry.build()
            except Exception as e:
                raise ValueError(f"Geometry '{name}' not built and auto-build failed: {e}. Call build() first.")

        transform = Transform3D(
            translation=position or (0.0, 0.0, 0.0),
            rotation=rotation or (0.0, 0.0, 0.0)
        )

        entry = ComponentEntry(
            geometry=geometry,
            key=key,
            base_name=base_name,
            transform=transform,
            aligned_port=align_port,
            metadata=metadata
        )

        self._components[key] = entry

        # Handle ordering - resolve component and optional gap
        gap = 0.0
        if isinstance(after, tuple):
            ref_after, gap = after
        else:
            ref_after = after

        if isinstance(before, tuple):
            ref_before, gap = before
        else:
            ref_before = before

        resolved_after = self._resolve_component_ref(ref_after) if ref_after else None
        resolved_before = self._resolve_component_ref(ref_before) if ref_before else None

        if resolved_after is not None:
            idx = self._component_order.index(resolved_after) + 1
            self._component_order.insert(idx, key)

            self._connections.append(Connection(
                from_key=resolved_after,
                to_key=key,
                from_port='port2',
                to_port=align_port or 'port1',
                gap=gap
            ))

        elif resolved_before is not None:
            idx = self._component_order.index(resolved_before)
            self._component_order.insert(idx, key)

            self._connections.append(Connection(
                from_key=key,
                to_key=resolved_before,
                from_port=align_port or 'port2',
                to_port='port1',
                gap=-gap  # Before means moving current component backwards relative to target
            ))
        else:
            self._component_order.append(key)

        self._layout_computed = False
        self._is_built = False
        self.invalidate_tag()
        self._record(
            'add', 
            name=name, 
            geometry_type=type(geometry).__name__,
            geometry_history=geometry.get_history(),
            position=position,
            rotation=rotation,
            after=after,
            before=before,
            align_port=align_port,
            **metadata
        )

        return self

    def _rename_first_instance(self, base_name: str) -> None:
        """Rename the first instance of a component when duplicates are added."""
        if base_name not in self._components:
            return

        # The first instance needs to be renamed to base_name_1
        old_key = base_name
        new_key = f"{base_name}_1"

        # Update component entry
        entry = self._components[old_key]
        entry.key = new_key

        # Move in dict
        self._components[new_key] = entry
        del self._components[old_key]

        # Update order list
        idx = self._component_order.index(old_key)
        self._component_order[idx] = new_key

        # Update connections
        for conn in self._connections:
            if conn.from_key == old_key:
                conn.from_key = new_key
            if conn.to_key == old_key:
                conn.to_key = new_key

        # Update base name groups
        if base_name in self._base_name_groups:
            self._base_name_groups[base_name] = [
                new_key if k == old_key else k 
                for k in self._base_name_groups[base_name]
            ]

        # Set count to 1 (next will be _2)
        self._base_name_counts[base_name] = 1

    def _resolve_component_ref(self, ref: str) -> Optional[str]:
        """
        Resolve a component reference (base name or full key) to actual key.

        If ref is a base name with multiple instances, returns the last one.
        """
        if ref is None:
            return None

        # Direct match
        if ref in self._components:
            return ref

        # Try as base name - return last instance
        if ref in self._base_name_groups and self._base_name_groups[ref]:
            return self._base_name_groups[ref][-1]

        raise ValueError(f"Component '{ref}' not found")

    def get_identical_components(self) -> Dict[str, List[str]]:
        """
        Get groups of components that share the same base name (geometry).

        These components are candidates for solver optimization - compute
        the solution once and reuse for all identical instances.

        Returns
        -------
        dict
            {base_name: [key1, key2, ...]} for groups with more than one member
        """
        return {
            name: keys for name, keys in self._base_name_groups.items()
            if len(keys) > 1
        }

    def get_components_by_base_name(self, base_name: str) -> List[str]:
        """Get all component keys with a given base name."""
        return self._base_name_groups.get(base_name, [])

    # =========================================================================
    # Connection and Layout Methods (unchanged logic, updated for new structure)
    # =========================================================================

    def connect(
        self,
        name: str,
        geometry: Union[BaseGeometry, 'Assembly'],
        to_component: str,
        from_port: str = "port1",
        to_port: str = "port2",
        rotation: Optional[Tuple[float, float, float]] = None,
        gap: float = 0.0,
        **metadata
    ) -> 'Assembly':
        """Add component connected via ports to existing component."""
        resolved_to = self._resolve_component_ref(to_component)

        self.add(name, geometry, rotation=rotation, **metadata)

        # Get the key that was just added (might be suffixed)
        new_key = self._component_order[-1]

        self._record(
            'connect',
            name=name,
            geometry_type=type(geometry).__name__,
            geometry_history=geometry.get_history(),
            to_component=to_component,
            from_port=from_port,
            to_port=to_port,
            rotation=rotation,
            gap=gap,
            **metadata
        )

        return self

    def remove(self, ref: str) -> 'Assembly':
        """Remove a component by key or base name (removes all instances if base name)."""
        keys_to_remove = []

        if ref in self._components:
            keys_to_remove = [ref]
        elif ref in self._base_name_groups:
            keys_to_remove = list(self._base_name_groups[ref])
        else:
            raise KeyError(f"Component '{ref}' not found")

        for key in keys_to_remove:
            entry = self._components[key]
            base_name = entry.base_name

            del self._components[key]
            self._component_order.remove(key)
            self._connections = [
                c for c in self._connections 
                if c.from_key != key and c.to_key != key
            ]

            # Update base name tracking
            if base_name in self._base_name_groups:
                self._base_name_groups[base_name].remove(key)
                if not self._base_name_groups[base_name]:
                    del self._base_name_groups[base_name]
                    del self._base_name_counts[base_name]

        self._layout_computed = False
        self._is_built = False
        self.invalidate_tag()

        self._record('remove', ref=ref)

        return self

    def rotate(
        self, 
        ref: str, 
        angles: Tuple[float, float, float],
        center: Optional[Tuple[float, float, float]] = None
    ) -> 'Assembly':
        """Apply rotation to a component."""
        key = self._resolve_component_ref(ref)
        entry = self._components[key]
        current = entry.transform.rotation
        entry.transform = Transform3D(
            translation=entry.transform.translation,
            rotation=tuple(c + a for c, a in zip(current, angles)),
            rotation_center=center
        )
        self._layout_computed = False
        self._is_built = False
        return self

    def translate(self, ref: str, offset: Tuple[float, float, float]) -> 'Assembly':
        """Apply translation to a component."""
        key = self._resolve_component_ref(ref)
        entry = self._components[key]
        current = entry.transform.translation
        entry.transform = Transform3D(
            translation=tuple(c + o for c, o in zip(current, offset)),
            rotation=entry.transform.rotation,
            rotation_center=entry.transform.rotation_center
        )
        self._layout_computed = False
        self._is_built = False
        return self

    # =========================================================================
    # Index-based Access
    # =========================================================================

    def __getitem__(self, key: Union[str, int]) -> ComponentEntry:
        """Access component by key, base name, or index."""
        if isinstance(key, int):
            if key < 0 or key >= len(self._component_order):
                raise IndexError(f"Index {key} out of range")
            key = self._component_order[key]
        else:
            key = self._resolve_component_ref(key)

        return self._components[key]

    def __len__(self) -> int:
        return len(self._components)

    def __iter__(self):
        return iter(self._component_order)

    def __contains__(self, ref: str) -> bool:
        try:
            self._resolve_component_ref(ref)
            return True
        except (KeyError, ValueError):
            return False

    @property
    def keys(self) -> List[str]:
        return list(self._component_order)

    @property
    def components(self) -> Dict[str, ComponentEntry]:
        return dict(self._components)

    # =========================================================================
    # Layout Computation (unchanged)
    # =========================================================================

    def compute_layout(self) -> 'Assembly':
        """Compute positions based on connections."""
        if self._layout_computed:
            return self

        if not self._components:
            self._layout_computed = True
            return self

        positioned = set()
        first_key = self._component_order[0]
        positioned.add(first_key)

        for conn in self._connections:
            if conn.to_key in positioned and conn.from_key in positioned:
                continue
            if conn.from_key in positioned:
                self._position_connected(conn)
                positioned.add(conn.to_key)
            elif conn.to_key in positioned:
                # Reverse connection positioning not yet implemented for all cases
                pass

        # Safety for unpositioned components
        for key in self._component_order:
            if key not in positioned:
                if positioned:
                    last = list(positioned)[-1]
                    self._position_after(key, last)
                positioned.add(key)

        self._layout_computed = True
        return self

    def _get_port_position(
        self, 
        entry: ComponentEntry, 
        port_name: str
    ) -> Optional[Tuple[float, float, float]]:
        """Get port position in local coordinates using physical bounds."""
        # Refresh physical bounds for better precision if possible
        if not isinstance(entry.geometry, Assembly) and entry.geometry.geo is not None:
            axis = self.main_axis
            z_min, z_max = entry.geometry.get_physical_bounds('Z')
            x_min, x_max = entry.geometry.get_physical_bounds('X')
            y_min, y_max = entry.geometry.get_physical_bounds('Y')
            pmin = (x_min, y_min, z_min)
            pmax = (x_max, y_max, z_max)
        else:
            if entry.original_bounds is None:
                return None
            pmin, pmax = entry.original_bounds

        axis_idx = self._AXIS_IDX[self.main_axis]
        center = [(pmin[i] + pmax[i]) / 2 for i in range(3)]

        if 'port1' in port_name.lower() or 'min' in port_name.lower():
            center[axis_idx] = pmin[axis_idx]
        elif 'port2' in port_name.lower() or 'max' in port_name.lower():
            center[axis_idx] = pmax[axis_idx]

        return tuple(center)

    def _position_connected(self, conn: Connection) -> None:
        """Position component based on port connection."""
        from_entry = self._components[conn.from_key]
        to_entry = self._components[conn.to_key]

        from_port_pos = self._get_port_position(from_entry, conn.from_port)
        to_port_pos = self._get_port_position(to_entry, conn.to_port)

        if from_port_pos is None or to_port_pos is None:
            self._position_after(conn.to_key, conn.from_key)
            return

        axis_idx = self._AXIS_IDX[self.main_axis]

        from_world = tuple(
            from_port_pos[i] + from_entry.transform.translation[i]
            for i in range(3)
        )

        translation = [0.0, 0.0, 0.0]

        for i in range(3):
            if i == axis_idx:
                translation[i] = from_world[i] - to_port_pos[i] + conn.gap
            else:
                from_center = (from_entry.original_bounds[0][i] + 
                             from_entry.original_bounds[1][i]) / 2
                to_center = (to_entry.original_bounds[0][i] + 
                           to_entry.original_bounds[1][i]) / 2
                translation[i] = (from_center + from_entry.transform.translation[i] 
                                 - to_center)

        to_entry.transform = Transform3D(
            translation=tuple(translation),
            rotation=to_entry.transform.rotation,
            rotation_center=to_entry.transform.rotation_center
        )

    def _position_after(self, key: str, after_key: str, gap: float = 0.0) -> None:
        """Position component after another along main axis using bounding boxes."""
        entry = self._components[key]
        after_entry = self._components[after_key]

        axis_idx = self._AXIS_IDX[self.main_axis]

        # World max of previous component
        pmin_after, pmax_after = after_entry.original_bounds
        world_max_after = pmax_after[axis_idx] + after_entry.transform.translation[axis_idx]

        # Local min of current component
        pmin_self, pmax_self = entry.original_bounds
        local_min_self = pmin_self[axis_idx]

        translation = [0.0, 0.0, 0.0]

        for i in range(3):
            if i == axis_idx:
                translation[i] = world_max_after - local_min_self + gap
            else:
                after_center = (pmin_after[i] + pmax_after[i]) / 2 + after_entry.transform.translation[i]
                self_center = (pmin_self[i] + pmax_self[i]) / 2
                translation[i] = after_center - self_center

        entry.transform = Transform3D(
            translation=tuple(translation),
            rotation=entry.transform.rotation,
            rotation_center=entry.transform.rotation_center
        )

    def _position_before(self, key: str, before_key: str, gap: float = 0.0) -> None:
        """Position component before another along main axis using bounding boxes."""
        entry = self._components[key]
        before_entry = self._components[before_key]

        axis_idx = self._AXIS_IDX[self.main_axis]

        # World min of next component
        pmin_before, pmax_before = before_entry.original_bounds
        world_min_before = pmin_before[axis_idx] + before_entry.transform.translation[axis_idx]

        # Local max of current component
        pmin_self, pmax_self = entry.original_bounds
        local_max_self = pmax_self[axis_idx]

        translation = [0.0, 0.0, 0.0]

        for i in range(3):
            if i == axis_idx:
                translation[i] = world_min_before - local_max_self - gap
            else:
                before_center = (pmin_before[i] + pmax_before[i]) / 2 + before_entry.transform.translation[i]
                self_center = (pmin_self[i] + pmax_self[i]) / 2
                translation[i] = before_center - self_center

        entry.transform = Transform3D(
            translation=tuple(translation),
            rotation=entry.transform.rotation,
            rotation_center=entry.transform.rotation_center
        )

    # =========================================================================
    # Build - Multi-solid Assembly
    # =========================================================================

    def build(self) -> None:
        """
        Build multi-solid geometry from all components.

        Creates a compound where each component is a separate solid with:
        - Solid named after component key (e.g., 'midcell_1', 'midcell_2')
        - Sequential port naming (port1, port2, ...) along main axis
        - Interface ports (where solids meet) shared between adjacent components
        - External ports at assembly ends
        - All other faces named 'wall'
        """
        if not self._components:
            raise ValueError("No components in assembly")

        self.compute_layout()

        # NEW: Snap interfaces to compensate for CAD inaccuracies
        self._snap_interfaces()

        shapes = []
        # Build sub-assemblies first
        for entry in self._components.values():
            if entry.is_assembly and not entry.geometry._is_built:
                entry.geometry.build()

        # Collect transformed shapes
        shapes = []
        shape_keys = []  # Track which key each shape belongs to

        for key in self._component_order:
            entry = self._components[key]
            shape = entry.geometry.geo
            transformed = self._apply_transform(shape, entry.transform)
            shapes.append(transformed)
            shape_keys.append(key)

        # To prevent merging identical/co-planar solids, we assign unique 
        # temporary materials and use OCCGeometry(shapes).shape which is
        # often more robust than Glue() for preserving solid identities.
        if len(shapes) > 1:
            for i, shape in enumerate(shapes):
                shape.mat(f"__temp_solid_{i}")

        # Create multi-solid compound
        if len(shapes) == 1:
            self.geo = shapes[0]
        else:
            # OCCGeometry(list) performs gluing but preserves solid domains better
            occ_geo = OCCGeometry(shapes)
            self.geo = occ_geo.shape

        # Name solids and setup boundaries
        self._name_solids(shape_keys)
        self._setup_boundaries()

        self._is_built = True
        self._record('build')

    def generate_mesh(self, maxh=None, curve_order: int = 3) -> Mesh:
        """
        Generate mesh with support for per-component refinement.
        """
        if self.geo is None:
            self.build()

        # Apply per-component maxh if available
        solids = list(self.geo.solids)
        if len(solids) == len(self._component_order):
            for i, key in enumerate(self._component_order):
                entry = self._components[key]
                # Try to get maxh from component geometry
                comp_maxh = getattr(entry.geometry, 'maxh', None)
                if comp_maxh:
                    solids[i].maxh = comp_maxh

        # Call base generate_mesh which will create OCCGeometry(self.geo)
        # and respect the per-solid maxh settings.
        return super().generate_mesh(maxh=maxh, curve_order=curve_order)

    def _name_solids(self, shape_keys: List[str]) -> None:
        """Name each solid in the geometry based on component keys.

        Each assembly component may contain multiple mesh solids (e.g. after
        PEC subtraction).  All solids belonging to the same component are
        given the same material name so the solver treats each component as
        a single domain.
        """
        self._solid_info = {}
        axis_idx = self._AXIS_IDX[self.main_axis]

        try:
            solids = list(self.geo.solids)
        except AttributeError:
            # Single solid
            key = shape_keys[0]
            entry = self._components[key]
            self.geo.mat(key)
            self._solid_info[key] = {
                'material': key,
                'base_name': entry.base_name,
                'position': 0,
                'bounds': None,
                'is_identical': False,
            }
            return

        # Compute bounding box for each component along the main axis
        component_ranges = []
        for key in shape_keys:
            entry = self._components[key]
            trf = entry.transform
            offset = trf.translation[axis_idx]
            # Get the component's own extent along main axis
            geo = entry.geometry.geo
            try:
                comp_solids = list(geo.solids)
            except AttributeError:
                comp_solids = [geo]
            lo = min(s.bounding_box[0][axis_idx] for s in comp_solids) + offset
            hi = max(s.bounding_box[1][axis_idx] for s in comp_solids) + offset
            component_ranges.append((lo, hi, key))

        # Sort components along main axis
        component_ranges.sort(key=lambda x: (x[0] + x[1]) / 2)

        # Assign each mesh solid to the component whose range contains it.
        # Solid names are prefixed with the component key so the solver can
        # map mesh materials back to their owning component while preserving
        # per-solid material distinctions for CoefficientFunction assembly.
        def get_solid_center(solid):
            bb = solid.bounding_box
            return (bb[0][axis_idx] + bb[1][axis_idx]) / 2

        # domain_materials: component_key -> [mesh_material_name, ...]
        self._domain_materials: Dict[str, List[str]] = {key: [] for _, _, key in component_ranges}

        for idx, solid in enumerate(solids):
            center = get_solid_center(solid)
            best_key = None
            best_dist = float('inf')
            for lo, hi, key in component_ranges:
                if lo - 0.01 <= center <= hi + 0.01:
                    best_key = key
                    break
                mid = (lo + hi) / 2
                dist = abs(center - mid)
                if dist < best_dist:
                    best_dist = dist
                    best_key = key

            # Prefix original solid name with component key for uniqueness
            orig_name = solid.name if solid.name else f"solid"
            # Strip temporary assembly naming
            if orig_name.startswith('__temp_solid_'):
                orig_name = f"solid"
            mat_name = f"{best_key}/{orig_name}"
            solid.mat(mat_name)
            self._domain_materials[best_key].append(mat_name)

        # Print domain material mapping
        print(f"\nAssembly domain-material mapping:")
        for key, mats in self._domain_materials.items():
            print(f"  {key}: {mats}")

        # Store component-level info
        for lo, hi, key in component_ranges:
            entry = self._components[key]
            self._solid_info[key] = {
                'material': key,
                'base_name': entry.base_name,
                'position': (lo + hi) / 2,
                'bounds': (lo, hi),
                'is_identical': len(self._base_name_groups.get(entry.base_name, [])) > 1,
                'mesh_materials': self._domain_materials[key],
            }

    def _snap_interfaces(self, tolerance: float = 1e-1) -> None: # 0.1 meter = 10cm tolerance (very aggressive)
        """
        Perform a bit-perfect alignment pass to snap faces together.
        """
        if not self._component_order or len(self._component_order) < 2:
            return

        axis_idx = self._AXIS_IDX[self.main_axis]

        # We follow the component order and snap each to its predecessor
        for i in range(1, len(self._component_order)):
            to_key = self._component_order[i]
            from_key = self._component_order[i-1]

            from_entry = self._components[from_key]
            to_entry = self._components[to_key]

            if from_entry.geometry.geo is None or to_entry.geometry.geo is None:
                continue

            # Get BIT-PRECISE physical extremes from extreme faces
            from_extremes = from_entry.geometry.get_extreme_faces(self.main_axis)
            to_extremes = to_entry.geometry.get_extreme_faces(self.main_axis)

            if 'max' not in from_extremes or 'min' not in to_extremes:
                # Fallback to physical bounds if no extreme face identified
                _, pmax_from = from_entry.geometry.get_physical_bounds(self.main_axis)
                pmin_to, _ = to_entry.geometry.get_physical_bounds(self.main_axis)
            else:
                pmax_from = from_extremes['max'][0]
                pmin_to = to_extremes['min'][0]

            world_max_from = pmax_from + from_entry.transform.translation[axis_idx]
            world_min_to = pmin_to + to_entry.transform.translation[axis_idx]

            gap = world_min_to - world_max_from

            # SNAP SAFETY: If the gap is massive, something is wrong with face detection
            # We don't want to fly components across the world.
            # Usually a gap should be < 10% of component size for snapping.
            comp_size = max(from_entry.size[axis_idx], to_entry.size[axis_idx])
            if abs(gap) > comp_size * 0.5:
                # print(f"SNAP REJECTED: Gap too large ({gap:.3e}) for {to_key}")
                continue

            # Snap them if they are reasonably close
            if abs(gap) < tolerance:
                trans = list(to_entry.transform.translation)
                trans[axis_idx] -= gap
                to_entry.transform.translation = tuple(trans)
                self._layout_computed = True

    def _name_solids_by_position(self, solids: list) -> None:
        """Fallback: name solids by their position along main axis."""
        axis_idx = self._AXIS_IDX[self.main_axis]

        def get_solid_position(solid):
            bb = solid.bounding_box
            return (bb[0][axis_idx] + bb[1][axis_idx]) / 2

        solids_sorted = sorted(solids, key=get_solid_position)

        for i, solid in enumerate(solids_sorted):
            name = f"cell_{i + 1}"
            solid.mat(name)
            bb = solid.bounding_box
            self._solid_info[name] = {
                'material': name,
                'base_name': name,
                'position': get_solid_position(solid),
                'bounds': (tuple(bb[0]), tuple(bb[1])),
                'is_identical': False
            }

    def _setup_boundaries(self) -> None:
        """
        Setup boundary conditions with unified sequential port naming.

        - All flat faces perpendicular to main axis are ports
        - Faces at the same position share the same port name (interfaces)
        - Port numbers increase along the positive main axis direction
        - Non-port faces are named 'wall'
        """
        if self.geo is None:
            return

        axis_idx = self._AXIS_IDX[self.main_axis]
        tolerance = 1e-3 # 1mm grouping tolerance for ports

        # Step 1: Name all faces 'wall' first
        self._name_all_faces_wall()

        # Step 2: Identify all port faces and their positions
        port_faces = []

        try:
            for face in self.geo.faces:
                try:
                    fbb = face.bounding_box
                    face_extent = fbb[1][axis_idx] - fbb[0][axis_idx]

                    if face_extent < tolerance:
                        face_pos = (fbb[0][axis_idx] + fbb[1][axis_idx]) / 2
                        port_faces.append((face_pos, face))
                except Exception:
                    pass
        except Exception as e:
            print(f"Warning: Error identifying faces: {e}")
            self.bc = 'default'
            self._bc_explicitly_set = True
            return

        if not port_faces:
            print("Warning: No port faces found")
            self.bc = 'default'
            self._bc_explicitly_set = True
            return

        # Step 3: Group faces by position
        port_faces.sort(key=lambda x: x[0])

        port_groups = []
        for pos, face in port_faces:
            matched = False
            for group in port_groups:
                if abs(group['position'] - pos) < tolerance:
                    group['faces'].append(face)
                    matched = True
                    break

            if not matched:
                port_groups.append({
                    'position': pos,
                    'faces': [face]
                })

        # Step 4: Sort groups by position and assign sequential port names
        port_groups.sort(key=lambda g: g['position'])

        self._port_info = {}

        for port_num, group in enumerate(port_groups, start=1):
            port_name = f'port{port_num}'

            is_external = len(group['faces']) == 1
            is_interface = len(group['faces']) > 1

            for face in group['faces']:
                face.name = port_name
                face.col = (1, 0, 0)

            # Find which components connect at this port
            connected_components = self._find_components_at_position(
                group['position'], tolerance
            )

            self._port_info[port_name] = {
                'position': group['position'],
                'num_faces': len(group['faces']),
                'type': 'external' if is_external else 'interface',
                'axis': self.main_axis,
                'connected_components': connected_components
            }

        # Summary
        n_external = sum(1 for p in self._port_info.values() if p['type'] == 'external')
        n_interface = sum(1 for p in self._port_info.values() if p['type'] == 'interface')

        print(f"Port naming complete:")
        print(f"  Total ports: {len(self._port_info)}")
        print(f"  External ports: {n_external}")
        print(f"  Interface ports: {n_interface}")

        self.bc = 'default'
        self._bc_explicitly_set = True

    def _name_all_faces_wall(self) -> None:
        """Name every face 'wall' as default."""
        if self.geo is None:
            return

        try:
            for solid in self.geo.solids:
                for face in solid.faces:
                    face.name = 'default'
        except AttributeError:
            try:
                for face in self.geo.faces:
                    face.name = 'default'
            except AttributeError:
                pass

    def get_material(self, domain_name: str) -> dict:
        """Get material properties for a mesh material, delegating to its component.

        Mesh materials for assemblies are prefixed with the component key,
        e.g. ``cell1/ceramic_1``.  This method strips the prefix and queries
        the owning component's geometry for the material properties.
        """
        defaults = {"eps_r": 1.0, "mu_r": 1.0, "sigma": 0.0, "tan_delta": 0.0}

        # Try to find the owning component via the domain_materials mapping
        dm = getattr(self, '_domain_materials', {})
        for comp_key, mat_names in dm.items():
            if domain_name in mat_names:
                entry = self._components[comp_key]
                # Strip component prefix to get the original solid name
                orig_name = domain_name[len(comp_key) + 1:] if domain_name.startswith(comp_key + '/') else domain_name
                if hasattr(entry.geometry, 'get_material'):
                    return entry.geometry.get_material(orig_name)
                return dict(defaults)

        # Fallback: try component key as domain_name
        if domain_name in self._components:
            return dict(defaults)

        return dict(defaults)

    def _find_components_at_position(
        self, 
        position: float, 
        tolerance: float
    ) -> List[Tuple[str, str]]:
        """Find which components have ports at a given position."""
        axis_idx = self._AXIS_IDX[self.main_axis]
        connected = []

        for key in self._component_order:
            entry = self._components[key]
            if entry.original_bounds is None:
                continue

            pmin, pmax = entry.original_bounds
            t = entry.transform.translation

            world_min = pmin[axis_idx] + t[axis_idx]
            world_max = pmax[axis_idx] + t[axis_idx]

            if abs(world_min - position) < tolerance:
                connected.append((key, 'port1'))
            if abs(world_max - position) < tolerance:
                connected.append((key, 'port2'))

        return connected

    def _apply_transform(self, shape, transform: Transform3D):
        """Apply transformation to OCC shape."""
        if transform.is_identity():
            return shape

        result = shape

        rx, ry, rz = transform.rotation
        if abs(rx) > 1e-12 or abs(ry) > 1e-12 or abs(rz) > 1e-12:
            if transform.rotation_center:
                center = transform.rotation_center
            else:
                try:
                    bb = shape.bounding_box
                    center = tuple((bb[0][i] + bb[1][i]) / 2 for i in range(3))
                except:
                    center = (0, 0, 0)

            if abs(rz) > 1e-12:
                result = result.Rotate(Axes(center, Z), rz)
            if abs(ry) > 1e-12:
                result = result.Rotate(Axes(center, Y), ry)
            if abs(rx) > 1e-12:
                result = result.Rotate(Axes(center, X), rx)

        tx, ty, tz = transform.translation
        if abs(tx) > 1e-12 or abs(ty) > 1e-12 or abs(tz) > 1e-12:
            result = result.Move((tx, ty, tz))

        return result

    # =========================================================================
    # Port and Solid Information
    # =========================================================================

    def get_port_info(self) -> Dict[str, Dict]:
        """Get information about all ports in the assembly."""
        return dict(self._port_info) if self._port_info else {}

    def get_solid_info(self) -> Dict[str, Dict]:
        """Get information about all solids in the assembly."""
        return dict(self._solid_info) if self._solid_info else {}

    def get_external_ports(self) -> List[str]:
        """Get list of external port names."""
        return [name for name, info in self._port_info.items() 
                if info['type'] == 'external']

    def get_interface_ports(self) -> List[str]:
        """Get list of interface port names."""
        return [name for name, info in self._port_info.items() 
                if info['type'] == 'interface']

    def print_port_info(self) -> None:
        """Print detailed port information."""
        if not self._port_info:
            print("No port information available. Call build() first.")
            return

        print("\n" + "=" * 70)
        print("ASSEMBLY PORT MAP")
        print("=" * 70)
        print(f"Main axis: {self.main_axis}")
        print(f"Total ports: {len(self._port_info)}")
        print("-" * 70)

        for port_name in sorted(self._port_info.keys(), 
                                key=lambda x: int(x.replace('port', ''))):
            info = self._port_info[port_name]
            port_type = info['type'].upper()
            pos = info['position']
            n_faces = info['num_faces']

            type_marker = "◀▶" if info['type'] == 'external' else "◀┃▶"

            components = info.get('connected_components', [])
            comp_str = ", ".join(f"{k}.{p}" for k, p in components) if components else "?"

            print(f"  {port_name:8s} │ {self.main_axis}={pos:+.6f} │ "
                  f"{port_type:10s} │ {n_faces} face(s) {type_marker}")
            print(f"           │ Components: {comp_str}")

        print("=" * 70)

        # Visual schematic
        self._print_schematic()

    def _print_schematic(self) -> None:
        """Print a visual schematic of the assembly."""
        ports_sorted = sorted(
            self._port_info.items(), 
            key=lambda x: x[1]['position']
        )

        print(f"\nSchematic (along {self.main_axis}):\n")

        # Build port line
        port_parts = []
        for port_name, info in ports_sorted:
            if info['type'] == 'external':
                port_parts.append(f"║{port_name}║")
            else:
                port_parts.append(f"┃{port_name}┃")

        # Build component line
        comp_parts = []
        for key in self._component_order:
            entry = self._components[key]
            if entry.base_name != key:
                # Show base name for duplicates
                comp_parts.append(f"[{key}]")
            else:
                comp_parts.append(f"[{key}]")

        print("  Ports:      " + " ─── ".join(port_parts))
        print("  Components: " + " ─── ".join(comp_parts))

        # Show identical component groups
        identical = self.get_identical_components()
        if identical:
            print(f"\n  Identical components (compute once, reuse):")
            for base_name, keys in identical.items():
                print(f"    {base_name}: {', '.join(keys)}")
        print()

    # =========================================================================
    # Visualization
    # =========================================================================

    def inspect(
        self,
        show_labels: bool = True,
        color_by_component: bool = True,
        **kwargs
    ) -> None:
        """Visualize assembly configuration."""
        self.compute_layout()

        if not self._components:
            print("Assembly is empty")
            return

        shapes = []
        colors = self._generate_colors(len(self._components))

        for i, key in enumerate(self._component_order):
            entry = self._components[key]
            shape = entry.geometry.geo
            transformed = self._apply_transform(shape, entry.transform)

            if color_by_component:
                try:
                    for face in transformed.faces:
                        face.col = colors[i]
                except:
                    pass

            shapes.append(transformed)

        if len(shapes) == 1:
            display_geo = shapes[0]
        else:
            display_geo = Glue(shapes)

        Draw(display_geo, **kwargs)

        # Summary
        print(f"\nAssembly: {len(self._components)} components")
        print("-" * 60)
        for i, key in enumerate(self._component_order):
            entry = self._components[key]
            pos = entry.transform.translation
            geo_type = type(entry.geometry).__name__
            size = entry.size

            base_info = f" (={entry.base_name})" if entry.base_name != key else ""
            sub = " [sub-assembly]" if entry.is_assembly else ""
            cached = " [CACHED]" if entry.geometry.has_cached_solution() else ""

            print(f"  [{i}] {key}{base_info}: {geo_type}{sub}{cached}")
            print(f"      Position: ({pos[0]:.4f}, {pos[1]:.4f}, {pos[2]:.4f})")
            print(f"      Size:     ({size[0]:.4f}, {size[1]:.4f}, {size[2]:.4f})")

        # Show identical groups
        identical = self.get_identical_components()
        if identical:
            print(f"\nIdentical components:")
            for base_name, keys in identical.items():
                print(f"  {base_name}: {len(keys)} instances -> compute once")

    def show(
        self,
        what: Literal["geometry", "mesh", "geo", "inspect"] = "geometry",
        **kwargs
    ) -> None:
        """Display the assembly."""
        what = what.lower()

        if what == "inspect":
            self.inspect(**kwargs)
            return
        elif what in ("geometry", "geo"):
            if self.geo is None:
                raise ValueError("Assembly not built. Call build() first.")
            scene = Draw(self.geo, **kwargs)
        elif what == "mesh":
            if self.mesh is None:
                raise ValueError("Mesh not generated. Call generate_mesh() first.")
            scene = Draw(self.mesh, **kwargs)
        else:
            raise ValueError(f"Invalid option '{what}'")

        _display_webgui_fallback(scene)

    @staticmethod
    def _generate_colors(n: int) -> List[Tuple[float, float, float]]:
        """Generate distinct colors."""
        colors = []
        for i in range(n):
            hue = i / max(n, 1)
            h = hue * 6
            c = 0.9 * 0.7
            x = c * (1 - abs(h % 2 - 1))
            m = 0.9 - c

            if h < 1: r, g, b = c, x, 0
            elif h < 2: r, g, b = x, c, 0
            elif h < 3: r, g, b = 0, c, x
            elif h < 4: r, g, b = 0, x, c
            elif h < 5: r, g, b = x, 0, c
            else: r, g, b = c, 0, x

            colors.append((r + m, g + m, b + m))
        return colors

    # =========================================================================
    # Assembly Bounds and Info
    # =========================================================================

    def get_assembly_bounds(self) -> Tuple[Tuple[float, ...], Tuple[float, ...]]:
        """Get bounding box of entire assembly."""
        if not self._components:
            return ((0, 0, 0), (0, 0, 0))

        all_min = [float('inf')] * 3
        all_max = [float('-inf')] * 3

        for key in self._components:
            entry = self._components[key]
            pmin, pmax = entry.original_bounds or ((0, 0, 0), (1, 1, 1))
            t = entry.transform.translation

            for i in range(3):
                all_min[i] = min(all_min[i], pmin[i] + t[i])
                all_max[i] = max(all_max[i], pmax[i] + t[i])

        return tuple(all_min), tuple(all_max)

    def summary(self) -> Dict[str, Any]:
        """Get assembly summary information."""
        identical = self.get_identical_components()

        return {
            'n_components': len(self._components),
            'n_unique_geometries': len(self._base_name_groups),
            'n_identical_groups': len(identical),
            'n_connections': len(self._connections),
            'n_ports': len(self._port_info),
            'n_external_ports': len(self.get_external_ports()),
            'n_interface_ports': len(self.get_interface_ports()),
            'main_axis': self.main_axis,
            'is_built': self._is_built,
            'has_mesh': self.mesh is not None,
            'component_keys': list(self._component_order),
            'identical_groups': identical,
            'bounds': self.get_assembly_bounds() if self._components else None,
        }

    def print_info(self) -> None:
        """Print detailed assembly information."""
        info = self.summary()

        print("\n" + "=" * 70)
        print("ASSEMBLY INFORMATION")
        print("=" * 70)
        print(f"Total components:       {info['n_components']}")
        print(f"Unique geometries:      {info['n_unique_geometries']}")
        print(f"Identical groups:       {info['n_identical_groups']}")
        print(f"Main axis:              {info['main_axis']}")
        print(f"Built:                  {info['is_built']}")
        print(f"Has mesh:               {info['has_mesh']}")

        print(f"\nPorts:")
        print(f"  Total:                {info['n_ports']}")
        print(f"  External:             {info['n_external_ports']}")
        print(f"  Interfaces:           {info['n_interface_ports']}")

        if info['bounds']:
            pmin, pmax = info['bounds']
            print(f"\nBounding Box:")
            print(f"  Min: ({pmin[0]:.4f}, {pmin[1]:.4f}, {pmin[2]:.4f})")
            print(f"  Max: ({pmax[0]:.4f}, {pmax[1]:.4f}, {pmax[2]:.4f})")

        print(f"\nComponents:")
        for i, key in enumerate(self._component_order):
            entry = self._components[key]
            geo_type = type(entry.geometry).__name__
            pos = entry.transform.translation

            base_info = f" (base: {entry.base_name})" if entry.base_name != key else ""

            print(f"  [{i}] {key}{base_info}: {geo_type}")
            print(f"       Position: ({pos[0]:.4f}, {pos[1]:.4f}, {pos[2]:.4f})")

        # Identical component groups
        if info['identical_groups']:
            print(f"\nIdentical Component Groups (solver optimization):")
            for base_name, keys in info['identical_groups'].items():
                print(f"  '{base_name}': {len(keys)} instances")
                print(f"    Keys: {', '.join(keys)}")
                print(f"    -> Compute once, reuse {len(keys)-1} times")

        print("=" * 70)

        # Print port map if built
        if self._is_built and self._port_info:
            self.print_port_info()

    # =========================================================================
    # Tagging and Cache Support
    # =========================================================================

    def _get_geometry_params(self) -> Dict[str, Any]:
        """Get parameters for tag computation."""
        component_tags = {}
        for key, entry in self._components.items():
            component_tags[key] = {
                'tag_hash': entry.tag.full_hash,
                'base_name': entry.base_name,
                'transform': {
                    'translation': entry.transform.translation,
                    'rotation': entry.transform.rotation
                }
            }

        return {
            'class': 'Assembly',
            'main_axis': self.main_axis,
            'components': component_tags,
            'order': list(self._component_order),
            'identical_groups': self.get_identical_components()
        }

    def _get_mesh_params(self) -> Dict[str, Any]:
        """Get mesh parameters for tag computation."""
        return {
            'maxh': getattr(self, 'maxh', None),
            'component_count': len(self._components)
        }

    def get_solver_optimization_info(self) -> Dict[str, Any]:
        """
        Get information for solver optimization.

        Returns details about which components are identical and can
        share computed solutions.

        Returns
        -------
        dict
            {
                'total_components': int,
                'unique_computations': int,
                'reusable_results': int,
                'groups': {
                    base_name: {
                        'keys': [key1, key2, ...],
                        'compute_key': key1,  # Which one to compute
                        'reuse_keys': [key2, ...]  # Which ones reuse the result
                    }
                }
            }
        """
        identical = self.get_identical_components()

        # Single-instance components
        unique_keys = [
            keys[0] for name, keys in self._base_name_groups.items()
            if len(keys) == 1
        ]

        groups = {}
        total_reuse = 0

        for base_name, keys in identical.items():
            groups[base_name] = {
                'keys': keys,
                'compute_key': keys[0],
                'reuse_keys': keys[1:]
            }
            total_reuse += len(keys) - 1

        return {
            'total_components': len(self._components),
            'unique_computations': len(self._base_name_groups),
            'reusable_results': total_reuse,
            'efficiency': 1 - (len(self._base_name_groups) / len(self._components)) 
                         if self._components else 0,
            'unique_keys': unique_keys,
            'groups': groups
        }

    def save_geometry(self, project_path) -> None:
        """
        Save assembly geometry and all sub-component source files.

        Creates:
        - ``geometry/components/{key}_source.{ext}`` for imported CAD files
        - ``geometry/components/{key}.step`` for primitives
        - ``geometry/history.json`` with rewritten paths
        """
        project_path = Path(project_path)
        geo_dir = project_path / 'geometry'
        comp_dir = geo_dir / 'components'
        geo_dir.mkdir(parents=True, exist_ok=True)
        comp_dir.mkdir(parents=True, exist_ok=True)

        # Also export the whole assembly as a STEP file for external tools
        if self.geo is not None:
            try:
                self.geo.WriteStep(str(geo_dir / 'assembly.step'))
            except Exception:
                pass

        # Build history with rewritten sub-component paths
        history_for_save = []
        component_sources = {}  # key -> {source_link, source_hash, source_filename}

        for entry in self._history:
            e = dict(entry)

            if e.get('op') == 'add' or e.get('op') == 'connect':
                name = e.get('name', '')
                geo_type = e.get('geometry_type', '')
                sub_history = e.get('geometry_history', [])

                # Find the component key for this name
                # (might be suffixed like name_1, name_2)
                comp_key = None
                for k, comp in self._components.items():
                    if comp.base_name == name:
                        if k not in component_sources:
                            comp_key = k
                            break
                if comp_key is None:
                    comp_key = name

                comp_entry = self._components.get(comp_key)
                if comp_entry is not None:
                    geo_obj = comp_entry.geometry
                    from .importers import OCCImporter

                    if isinstance(geo_obj, OCCImporter) and geo_obj._source_link:
                        # Copy source CAD file
                        src = Path(geo_obj._source_link)
                        if src.exists():
                            ext = src.suffix
                            local_name = f'{comp_key}_source{ext}'
                            dest = comp_dir / local_name
                            if src.absolute() != dest.absolute():
                                shutil.copy2(str(src), str(dest))
                            source_hash = self._file_hash(dest)
                            component_sources[comp_key] = {
                                'source_link': str(geo_obj._source_link),
                                'source_hash': source_hash,
                                'source_filename': f'components/{local_name}',
                            }
                            # Rewrite filepath in sub-history
                            rewritten_sub = []
                            for sh in sub_history:
                                sh_copy = dict(sh)
                                if sh_copy.get('op') in ('import_occ', 'import_step'):
                                    sh_copy['filepath'] = f'geometry/components/{local_name}'
                                rewritten_sub.append(sh_copy)
                            e['geometry_history'] = rewritten_sub

                    elif not isinstance(geo_obj, Assembly):
                        # Primitive — export as STEP
                        if geo_obj.geo is not None:
                            try:
                                step_name = f'{comp_key}.step'
                                geo_obj.geo.WriteStep(str(comp_dir / step_name))
                                component_sources[comp_key] = {
                                    'source_link': None,
                                    'source_hash': None,
                                    'source_filename': f'components/{step_name}',
                                }
                            except Exception:
                                pass

            history_for_save.append(e)

        meta = {
            'type': 'Assembly',
            'module': 'cavsim3d.geometry.assembly',
            'source_link': None,
            'source_filename': None,
            'source_hash': None,
            'component_sources': component_sources,
            'history': history_for_save,
        }

        with open(geo_dir / 'history.json', 'w') as f:
            json.dump(meta, f, indent=2, default=str)

    @classmethod
    def _rebuild_from_history(
        cls,
        history: List[dict],
        project_path,
        source_file=None,
    ) -> 'Assembly':
        """Reconstruct assembly by replaying operation history."""
        project_path = Path(project_path)
        obj = None

        for entry in history:
            op = entry['op']
            params = {k: v for k, v in entry.items() if k not in ['op', 'timestamp']}

            if op == '__init__':
                obj = cls(main_axis=params.get('main_axis', 'Z'))

            elif op == 'add':
                # Reconstruct sub-geometry
                sub_history = list(params.pop('geometry_history'))
                sub_type = params.pop('geometry_type')

                sub_cls = BaseGeometry._get_subclass(sub_type)
                if sub_cls is None:
                    raise ValueError(f"Unknown geometry type '{sub_type}'")

                # Resolve source file from rewritten paths
                sub_source = None
                for sh in sub_history:
                    if sh.get('op') in ('import_occ', 'import_step'):
                        fp = sh.get('filepath', '')
                        candidate = project_path / fp
                        if candidate.exists():
                            sub_source = candidate
                        break

                sub_geo = sub_cls._rebuild_from_history(
                    sub_history, project_path, source_file=sub_source
                )

                name = params.pop('name')
                obj.add(name, sub_geo, **params)

            elif op == 'connect':
                sub_history = list(params.pop('geometry_history'))
                sub_type = params.pop('geometry_type')

                sub_cls = BaseGeometry._get_subclass(sub_type)
                if sub_cls is None:
                    raise ValueError(f"Unknown geometry type '{sub_type}'")

                sub_source = None
                for sh in sub_history:
                    if sh.get('op') in ('import_occ', 'import_step'):
                        fp = sh.get('filepath', '')
                        candidate = project_path / fp
                        if candidate.exists():
                            sub_source = candidate
                        break

                sub_geo = sub_cls._rebuild_from_history(
                    sub_history, project_path, source_file=sub_source
                )

                name = params.pop('name')
                obj.connect(name, sub_geo, **params)

            elif op == 'remove':
                obj.remove(**params)

            elif op == 'rotate':
                obj.rotate(**params)

            elif op == 'translate':
                obj.translate(**params)

            elif op == 'build':
                obj.build()

            elif op == 'generate_mesh':
                obj.generate_mesh(**params)

        return obj

components property