Coverage for fiqus/geom_generators/GeometryConductorAC_CC.py: 82%
825 statements
« prev ^ index » next coverage.py v7.4.4, created at 2026-01-09 01:49 +0000
1"""
2GeometryConductorAC_CC.py
42D rectangular stack geometry for a striated coated-conductor model.
6Strategy
7-------------------
8- Build a stack of rectangular layers aligned along x (width) and y (thickness):
9 HTS
10 Substrate
11 optional Silver (top / bottom)
12 optional Copper (top / bottom / left / right)
13- Each layer is created as a Gmsh OpenCASCADE rectangle (surface).
14- For the HTS layer we optionally striate the superconductor into
15 multiple filaments separated by grooves, using rectangle cuts.
16- After all layers are built:
17 - We create a circular air region around the stack.
18 - We classify boundary edges of each layer into Upper / Lower /
19 LeftEdge / RightEdge using only bounding boxes (robust to OCC ops).
20 - We create a consistent set of 2D/1D physical groups for
21 FiQuS/GetDP (materials + interfaces + outer air boundary).
23surface_tag = OCC surface entity (dim=2)
24curve_tag = boundary curve entity (dim=1)
25"""
27import logging
28import math
29import os
30from types import SimpleNamespace
31from typing import Dict, Iterable, List, Optional, Tuple
33import gmsh
35import fiqus.data.DataFiQuSConductor as geom # TODO
36from fiqus.data.DataFiQuS import FDM
37from fiqus.utils.Utils import GmshUtils, FilesAndFolders
38from fiqus.data.DataConductor import CC, Conductor
40logger = logging.getLogger("FiQuS")
42# ---------------------------------------------------------------------------
43# Small validation helpers
44# ---------------------------------------------------------------------------
46def _as_float(value, *, name: str) -> float:
47 try:
48 out = float(value)
49 except (TypeError, ValueError) as exc:
50 raise ValueError(f"{name} must be a real number, got {value!r}") from exc
51 if not math.isfinite(out):
52 raise ValueError(f"{name} must be finite, got {out!r}")
53 return out
56def _require_positive(value, *, name: str) -> float:
57 v = _as_float(value, name=name)
58 if v <= 0.0:
59 raise ValueError(f"{name} must be > 0, got {v:g}")
60 return v
63def _require_non_negative_optional(value, *, name: str) -> float:
64 """
65 Optional thickness rule:
66 - None => treated as 0 (layer absent)
67 - >= 0 => ok
68 - < 0 => hard error
69 """
70 if value is None:
71 return 0.0
72 v = _as_float(value, name=name)
73 if v < 0.0:
74 raise ValueError(f"{name} must be >= 0, got {v:g}")
75 return v
78def _require_int_ge(value, *, name: str, min_value: int) -> Optional[int]:
79 if value is None:
80 return None
81 try:
82 iv = int(value)
83 except (TypeError, ValueError) as exc:
84 raise ValueError(f"{name} must be an integer, got {value!r}") from exc
85 if iv < min_value:
86 raise ValueError(f"{name} must be >= {min_value}, got {iv}")
87 return iv
90# ---------------------------------------------------------------------------
91# Base class: axis-aligned rectangular layer
92# ---------------------------------------------------------------------------
94class _RectLayerBase:
95 """
96 Common base for all axis-aligned rectangular layers (HTS, substrate,
97 silver, copper).
99 Responsibilities
100 ----------------
101 - Store the main OCC surface (surface_tag).
102 - Maintain lists of boundary curves (curve_tags) and points (point_tags).
103 - Classify boundary curves into semantic edges:
104 "Upper", "Lower", "LeftEdge", "RightEdge"
105 based solely on bounding boxes (robust after OCC boolean ops).
106 """
108 def __init__(self) -> None:
109 self.surface_tag: Optional[int] = None
110 self.curve_tags: List[int] = []
111 self.point_tags: List[int] = []
112 self.edge_tags: Dict[str, List[int]] = {}
114 # Basic OCC -> topology refresh
116 def _refresh_from_surface(self) -> None:
117 """
118 Refresh curve_tags and point_tags from the current surface_tag.
119 This should be called after any OCC operation (cut, fuse, etc.)
120 that may re-tag curves/points.
121 """
122 if self.surface_tag is None:
123 raise RuntimeError(
124 f"{self.__class__.__name__} has no surface_tag to refresh from."
125 )
127 boundary = gmsh.model.getBoundary( [(2, self.surface_tag)], oriented=False, recursive=False )
128 self.curve_tags = [c[1] for c in boundary]
130 pts: List[int] = []
131 for c in self.curve_tags:
132 ends = gmsh.model.getBoundary( [(1, c)], oriented=False, recursive=False )
133 pts.extend([p[1] for p in ends])
135 # De-duplicate while preserving order
136 self.point_tags = list(dict.fromkeys(pts))
138 # Edge classification
140 def _classify_edges_from_bbox(self, w: float, t: float, cx: float, cy: float ) -> None:
141 """
142 Classify boundary curves into:
143 Upper / Lower / LeftEdge / RightEdge
145 Strategy
146 --------
147 - For each boundary curve, read its bounding box.
148 - Decide whether the curve is "horizontal" (dx >= dy) or
149 "vertical" (dy > dx).
150 - Among all horizontal curves:
151 - Those with y close to the maximum are "Upper".
152 - Those with y close to the minimum are "Lower".
153 - Among all vertical curves:
154 - Those with x close to the minimum are "LeftEdge".
155 - Those with x close to the maximum are "RightEdge".
156 - Tolerance is 5% of the global span in x/y for that set.
158 This is robust against most OCC operations as long as the
159 rectangle is still axis-aligned overall.
160 OCC fragment/cut operations can split edges into multiple curves and reorder
161 tags. Using curve bounding boxes lets us recover "Upper/Lower/Left/Right" as
162 long as the overall layer remains axis-aligned.
163 """
164 horiz: List[Tuple[int, float]] = [] # (curve_tag, y_mid)
165 vert: List[Tuple[int, float]] = [] # (curve_tag, x_mid)
167 for c in self.curve_tags:
168 xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(1, c)
169 dx = abs(xmax - xmin)
170 dy = abs(ymax - ymin)
172 if dx >= dy:
173 # Mostly horizontal
174 ymid = 0.5 * (ymin + ymax)
175 horiz.append((c, ymid))
176 else:
177 # Mostly vertical
178 xmid = 0.5 * (xmin + xmax)
179 vert.append((c, xmid))
181 if not horiz or not vert:
182 raise RuntimeError(
183 f"Could not find both horizontal and vertical edges for "
184 f"surface {self.surface_tag}; curves={self.curve_tags}"
185 )
187 # Horizontal: Upper / Lower
188 ys = [ym for _, ym in horiz]
189 y_min = min(ys)
190 y_max = max(ys)
191 span_y = max(y_max - y_min, 1e-12) # avoid zero span
192 tol_y = 0.05 * span_y # 5% band near extremes
194 upper = [c for c, ym in horiz if (y_max - ym) <= tol_y]
195 lower = [c for c, ym in horiz if (ym - y_min) <= tol_y]
197 # Vertical: Left / Right
198 xs = [xm for _, xm in vert]
199 x_min = min(xs)
200 x_max = max(xs)
201 span_x = max(x_max - x_min, 1e-12)
202 tol_x = 0.05 * span_x
204 left = [c for c, xm in vert if (xm - x_min) <= tol_x]
205 right = [c for c, xm in vert if (x_max - xm) <= tol_x]
207 edge_tags: Dict[str, List[int]] = {}
208 if upper:
209 edge_tags["Upper"] = upper
210 if lower:
211 edge_tags["Lower"] = lower
212 if left:
213 edge_tags["LeftEdge"] = left
214 if right:
215 edge_tags["RightEdge"] = right
217 missing = {"Upper", "Lower", "LeftEdge", "RightEdge"} - set(edge_tags)
218 if missing:
219 raise RuntimeError(
220 f"Edges not fully classified: {missing}; "
221 f"curves={self.curve_tags}, horiz={horiz}, vert={vert}"
222 )
224 self.edge_tags = edge_tags
226 # ---- public refresh API -----------------------------------------------
228 def refresh_topology(self, w: float, t: float, cx: float, cy: float) -> None:
229 """
230 Re-pull curves/points from the current surface and re-classify edges.
232 Subclasses that use multiple surfaces (e.g. striated HTS) can
233 override this method but should call _refresh_from_surface and
234 _classify_edges_from_bbox at some point.
235 """
236 self._refresh_from_surface()
237 self._classify_edges_from_bbox(w, t, cx, cy)
240# ---------------------------------------------------------------------------
241# HTS layer: optional striation into filaments + grooves
242# ---------------------------------------------------------------------------
244class HTSLayer(_RectLayerBase):
245 """
246 2D HTS layer, optionally striated into multiple filaments.
248 Geometric parameters
249 --------------------
250 - HTS_thickness, HTS_width
251 - HTS_center_x, HTS_center_y
253 Striation model
254 ---------------
255 If (n_striations > 1 and striation_w > 0):
256 - We build a full-width HTS rectangle.
257 - We build (N - 1) vertical groove rectangles.
258 - We OCC-cut the grooves out of the HTS.
259 - The remaining islands are HTS filaments.
260 - Grooves are kept as separate surfaces.
262 Bookkeeping
263 -----------
264 - self.surface_tag:
265 - monolithic HTS: the only HTS surface
266 - striated HTS: representative filament tag (for legacy code)
267 - self.filament_tags: list of all HTS filament surface IDs
268 - self.groove_tags: list of all groove surface IDs
269 """
271 def __init__(
272 self,
273 HTS_thickness: float,
274 HTS_width: float,
275 HTS_center_x: float,
276 HTS_center_y: float,
277 number_of_filaments: Optional[int] = None,
278 gap_between_filaments: Optional[float] = None,
279 ) -> None:
280 super().__init__()
282 self.HTS_thickness = float(HTS_thickness)
283 self.HTS_width = float(HTS_width)
284 self.HTS_center_x = float(HTS_center_x)
285 self.HTS_center_y = float(HTS_center_y)
287 # Striation parameters
288 self.n_striations = (
289 int(number_of_filaments) if number_of_filaments is not None else None
290 )
291 self.striation_w = (
292 float(gap_between_filaments)
293 if gap_between_filaments is not None
294 else None
295 )
297 # After cutting
298 self.filament_tags: List[int] = [] # HTS islands
299 self.groove_tags: List[int] = [] # Groove surfaces
304 def _refresh_from_surface(self) -> None:
305 """
306 Refresh boundary curves/points for the HTS layer.
308 Notes on striated HTS
309 ---------------------
310 When striation is enabled, the HTS is represented by multiple OCC surfaces
311 (one per filament). For boundary and interface physical groups we treat the
312 HTS as the *union of filaments* and recover its boundary by collecting the
313 boundaries of all filament surfaces.
315 This deliberately ignores groove surfaces: grooves are separate regions and
316 must not be part of the HTS conductor boundary.
318 Fallback
319 --------
320 If HTS is monolithic (no filaments), we fall back to the base implementation
321 using self.surface_tag.
322 """
324 if self.filament_tags:
325 curves: List[int] = []
327 for s in self.filament_tags:
328 boundary = gmsh.model.getBoundary( [(2, s)], oriented=False, recursive=False )
329 curves.extend( [c[1] for c in boundary] )
331 self.curve_tags = list(dict.fromkeys(curves))
333 pts: List[int] = []
334 for c in self.curve_tags:
335 ends = gmsh.model.getBoundary( [(1, c)], oriented=False, recursive=False )
336 pts.extend( [p[1] for p in ends] )
337 self.point_tags = list(dict.fromkeys(pts))
338 else:
339 super()._refresh_from_surface()
343 def build_HTS(self) -> int:
344 """
345 Build a 2D rectangular HTS layer centered at (HTS_center_x, HTS_center_y).
347 If striation parameters are valid, cut the HTS into multiple
348 filaments by removing narrow vertical grooves.
350 Returns
351 -------
352 int
353 A representative HTS surface tag (for backward compatibility).
354 """
355 x0 = self.HTS_center_x - self.HTS_width / 2.0
356 y0 = self.HTS_center_y - self.HTS_thickness / 2.0
358 # Base full-width HTS rectangle
359 base_tag = gmsh.model.occ.addRectangle( x0, y0, 0.0, self.HTS_width, self.HTS_thickness )
360 self.surface_tag = base_tag # initial monolithic tag
362 grooves_occ: List[Tuple[int, int]] = [] # [(2, tag), ...]
363 self.groove_tags = []
364 self.filament_tags = []
366 # Decide whether to striate
367 if (
368 self.n_striations is not None
369 and self.n_striations > 1
370 and self.striation_w is not None
371 and self.striation_w > 0.0
372 ):
373 N = self.n_striations
374 gw = self.striation_w
376 total_groove = gw * (N - 1)
377 if total_groove >= self.HTS_width:
378 logger.warning(
379 "[Geom-CC] Requested HTS striations do not fit in HTS width; "
380 "skipping striation (N=%d, gw=%.3g, W=%.3g).",
381 N,
382 gw,
383 self.HTS_width,
384 )
385 else:
386 # Ensure:
387 # N * filament_width + (N - 1) * gw = HTS_width
388 filament_w = (self.HTS_width - total_groove) / N
390 x_left = x0
391 for i in range(N - 1):
392 # Position of groove i between filament i and i+1
393 xg = x_left + (i + 1) * filament_w + i * gw
394 gtag = gmsh.model.occ.addRectangle( xg, y0, 0.0, gw, self.HTS_thickness )
395 grooves_occ.append((2, gtag))
396 self.groove_tags.append(gtag)
398 # Cut grooves out of HTS:
399 # removeObject=True -> base HTS replaced by its parts
400 # removeTool=False -> keep the groove rectangles as surfaces
401 cut_objs, _ = gmsh.model.occ.cut( [(2, base_tag)], grooves_occ, removeObject=True, removeTool=False )
403 if not cut_objs:
404 raise RuntimeError("[Geom-CC] HTS striation cut produced no surfaces.")
406 self.filament_tags = [e[1] for e in cut_objs]
407 # Representative tag (for legacy code paths)
408 self.surface_tag = self.filament_tags[0]
410 gmsh.model.occ.synchronize()
411 self.refresh_topology()
413 return self.surface_tag
416 def refresh_topology(self) -> None:
417 """
418 Refresh HTS topology and edge classification using global
419 HTS bounding box parameters.
420 """
421 super().refresh_topology(
422 self.HTS_width,
423 self.HTS_thickness,
424 self.HTS_center_x,
425 self.HTS_center_y,
426 )
428 # Physical groups
430 def create_physical_groups(self, name_prefix: str = "HTS") -> Dict[str, int]:
431 """
432 Create physical groups for the HTS layer, its filaments, grooves,
433 and edges.
435 2D groups
436 ---------
437 - name_prefix:
438 all HTS filaments together (monolithic or striated).
439 - name_prefix_Filament_i:
440 individual HTS filaments (ordered left-to-right).
441 - name_prefix_Grooves:
442 all grooves together (if any).
443 - name_prefix_Groove_i:
444 individual groove surfaces (ordered left-to-right).
446 1D edge groups
447 --------------
448 - name_prefix_Upper
449 - name_prefix_Lower
450 - name_prefix_LeftEdge
451 - name_prefix_RightEdge
452 """
453 if self.surface_tag is None and not self.filament_tags:
454 raise RuntimeError("build_HTS() must be called before create_physical_groups().")
456 out: Dict[str, int] = {}
458 # Helper for left-to-right ordering
459 def x_center_surf(tag: int) -> float:
460 xmin, _, _, xmax, _, _ = gmsh.model.getBoundingBox(2, tag)
461 return 0.5 * (xmin + xmax)
463 # Surfaces: all filaments together (striation-aware)
464 if self.filament_tags:
465 surf_tags = self.filament_tags
466 else:
467 surf_tags = [self.surface_tag]
469 pg = gmsh.model.addPhysicalGroup(2, surf_tags)
470 gmsh.model.setPhysicalName(2, pg, name_prefix)
471 out[name_prefix] = pg
473 # Per-filament surface groups
474 if self.filament_tags:
475 for idx, s in enumerate( sorted(self.filament_tags, key=x_center_surf), start=1 ):
476 pg_f = gmsh.model.addPhysicalGroup(2, [s])
477 name_f = f"{name_prefix}_Filament_{idx}"
478 gmsh.model.setPhysicalName(2, pg_f, name_f)
479 out[name_f] = pg_f
481 # Groove surface groups
482 if self.groove_tags:
483 # All grooves together
484 pg_all_g = gmsh.model.addPhysicalGroup(2, self.groove_tags)
485 name_all_g = f"{name_prefix}_Grooves"
486 gmsh.model.setPhysicalName(2, pg_all_g, name_all_g)
487 out[name_all_g] = pg_all_g
489 # Individual grooves: left -> right
490 for idx, s in enumerate( sorted(self.groove_tags, key=x_center_surf), start=1 ):
491 pg_g = gmsh.model.addPhysicalGroup(2, [s])
492 name_g = f"{name_prefix}_Groove_{idx}"
493 gmsh.model.setPhysicalName(2, pg_g, name_g)
494 out[name_g] = pg_g
496 # Edge groups (union over all filaments)
497 for name, tags in self.edge_tags.items():
498 edge_pg = gmsh.model.addPhysicalGroup(1, tags)
499 gmsh.model.setPhysicalName(1, edge_pg, f"{name_prefix}_{name}")
500 out[f"{name_prefix}_{name}"] = edge_pg
502 return out
505# ---------------------------------------------------------------------------
506# Substrate layer: underneath HTS
507# ---------------------------------------------------------------------------
509class SubstrateLayer(_RectLayerBase):
510 """
511 2D rectangular substrate, placed directly underneath the HTS layer.
513 - Shares a flat interface with the HTS bottom.
514 - Same width and x-center as HTS.
515 """
517 def __init__(self, substrate_thickness: float) -> None:
518 super().__init__()
520 self.substrate_thickness = float(substrate_thickness)
521 self.width: Optional[float] = None
522 self.center_x: Optional[float] = None
523 self.center_y: Optional[float] = None
525 def build_substrate(self, hts: HTSLayer) -> int:
526 """
527 Place the substrate directly under the provided HTS layer
528 (touching at the HTS bottom).
530 Parameters
531 ----------
532 hts : HTSLayer
533 The HTS layer that must already be built.
535 Returns
536 -------
537 int
538 The OCC surface tag of the substrate.
539 """
540 if hts.surface_tag is None:
541 raise RuntimeError(
542 "HTS layer must be built before building the substrate."
543 )
545 w = hts.HTS_width
546 t = hts.HTS_thickness
547 cx = hts.HTS_center_x
548 cy = hts.HTS_center_y - t / 2.0 - self.substrate_thickness / 2.0
550 self.width = w
551 self.center_x = cx
552 self.center_y = cy
554 x0 = cx - w / 2.0
555 y0 = cy - self.substrate_thickness / 2.0
557 self.surface_tag = gmsh.model.occ.addRectangle(x0, y0, 0.0, w, self.substrate_thickness)
559 gmsh.model.occ.synchronize()
560 self.refresh_topology()
561 return self.surface_tag
563 def refresh_topology(self) -> None:
564 if self.width is None or self.center_x is None or self.center_y is None:
565 raise RuntimeError("Substrate geometry not yet initialized.")
566 super().refresh_topology(
567 self.width,
568 self.substrate_thickness,
569 self.center_x,
570 self.center_y,
571 )
573 def create_physical_groups(self, name_prefix: str = "Substrate") -> Dict[str, int]:
574 """
575 Create physical groups for the substrate layer and its edges.
576 """
577 if self.surface_tag is None:
578 raise RuntimeError("build_substrate() must be called before create_physical_groups().")
580 out: Dict[str, int] = {}
582 pg = gmsh.model.addPhysicalGroup(2, [self.surface_tag])
583 gmsh.model.setPhysicalName(2, pg, name_prefix)
584 out[name_prefix] = pg
586 for name, tags in self.edge_tags.items():
587 edge_pg = gmsh.model.addPhysicalGroup(1, tags)
588 gmsh.model.setPhysicalName(1, edge_pg, f"{name_prefix}_{name}")
589 out[f"{name_prefix}_{name}"] = edge_pg
591 return out
594# ---------------------------------------------------------------------------
595# Silver layers: thin stabilizer above HTS / below substrate
596# ---------------------------------------------------------------------------
598class SilverTopLayer(_RectLayerBase):
599 """
600 Top silver layer placed directly above HTS (and directly below CopperTop).
602 This layer:
603 - Has same width and x-center as HTS.
604 - Touches the HTS upper surface (or can be removed via config).
605 """
607 def __init__(self, thickness: float) -> None:
608 super().__init__()
609 self.thickness = float(thickness)
610 self.width: Optional[float] = None
611 self.center_x: Optional[float] = None
612 self.center_y: Optional[float] = None
615 def build_over(self, hts: HTSLayer) -> int:
616 """
617 Place the top silver layer directly above the HTS layer.
618 """
619 if hts.surface_tag is None:
620 raise RuntimeError("HTS must be built before the top Silver layer.")
622 w = hts.HTS_width
623 cx = hts.HTS_center_x
624 cy = hts.HTS_center_y + hts.HTS_thickness / 2.0 + self.thickness / 2.0
626 self.width = w
627 self.center_x = cx
628 self.center_y = cy
630 x0 = cx - w / 2.0
631 y0 = cy - self.thickness / 2.0
633 self.surface_tag = gmsh.model.occ.addRectangle( x0, y0, 0.0, w, self.thickness )
635 gmsh.model.occ.synchronize()
636 self.refresh_topology()
637 return self.surface_tag
639 def refresh_topology(self) -> None:
640 if self.width is None or self.center_x is None or self.center_y is None:
641 raise RuntimeError("SilverTop geometry not yet initialized.")
642 super().refresh_topology(
643 self.width,
644 self.thickness,
645 self.center_x,
646 self.center_y,
647 )
649 def create_physical_groups(self, name_prefix: str = "SilverTop") -> Dict[str, int]:
650 """
651 Create physical groups for the top silver layer and its edges.
652 """
653 if self.surface_tag is None:
654 raise RuntimeError("build_over() must be called before creating PGs.")
656 out: Dict[str, int] = {}
658 pg = gmsh.model.addPhysicalGroup(2, [self.surface_tag])
659 gmsh.model.setPhysicalName(2, pg, name_prefix)
660 out[name_prefix] = pg
662 for name, tags in self.edge_tags.items():
663 epg = gmsh.model.addPhysicalGroup(1, tags)
664 gmsh.model.setPhysicalName(1, epg, f"{name_prefix}_{name}")
665 out[f"{name_prefix}_{name}"] = epg
667 return out
670class SilverBottomLayer(_RectLayerBase):
671 """
672 Bottom silver layer placed directly under the substrate
673 (and directly above CopperBottom).
674 """
676 def __init__(self, thickness: float) -> None:
677 super().__init__()
678 self.thickness = float(thickness)
679 self.width: Optional[float] = None
680 self.center_x: Optional[float] = None
681 self.center_y: Optional[float] = None
684 def build_under(self, sub: SubstrateLayer) -> int:
685 """
686 Place SilverBottom directly under the Substrate layer (touching
687 at the interface).
688 """
689 if sub.surface_tag is None:
690 raise RuntimeError("Substrate must be built before the bottom Silver layer.")
692 w = sub.width
693 cx = sub.center_x
694 cy = sub.center_y - sub.substrate_thickness / 2.0 - self.thickness / 2.0
696 if w is None or cx is None or cy is None:
697 raise RuntimeError("Substrate geometry not yet initialized.")
699 self.width = w
700 self.center_x = cx
701 self.center_y = cy
703 x0 = cx - w / 2.0
704 y0 = cy - self.thickness / 2.0
706 self.surface_tag = gmsh.model.occ.addRectangle(x0, y0, 0.0, w, self.thickness)
708 gmsh.model.occ.synchronize()
709 self.refresh_topology()
710 return self.surface_tag
712 def refresh_topology(self) -> None:
713 if self.width is None or self.center_x is None or self.center_y is None:
714 raise RuntimeError("SilverBottom geometry not yet initialized.")
715 super().refresh_topology(
716 self.width,
717 self.thickness,
718 self.center_x,
719 self.center_y,
720 )
722 def create_physical_groups(self, name_prefix: str = "SilverBottom") -> Dict[str, int]:
723 """
724 Create physical groups for the bottom silver layer and its edges.
725 """
726 if self.surface_tag is None:
727 raise RuntimeError("build_under() must be called before creating PGs.")
729 out: Dict[str, int] = {}
731 pg = gmsh.model.addPhysicalGroup(2, [self.surface_tag])
732 gmsh.model.setPhysicalName(2, pg, name_prefix)
733 out[name_prefix] = pg
735 for name, tags in self.edge_tags.items():
736 epg = gmsh.model.addPhysicalGroup(1, tags)
737 gmsh.model.setPhysicalName(1, epg, f"{name_prefix}_{name}")
738 out[f"{name_prefix}_{name}"] = epg
740 return out
743# ---------------------------------------------------------------------------
744# Copper layers: bottom / top / left / right
745# ---------------------------------------------------------------------------
747class CopperBottomLayer(_RectLayerBase):
748 """
749 Lower 2D copper layer placed directly under a base layer:
751 - Under SilverBottom if present;
752 - Otherwise directly under Substrate.
754 This is handled in the geometry builder (Generate_geometry).
755 """
757 def __init__(self, thickness: float) -> None:
758 super().__init__()
759 self.thickness = float(thickness)
760 self.width: Optional[float] = None
761 self.center_x: Optional[float] = None
762 self.center_y: Optional[float] = None
764 def build_under(self, base_layer: _RectLayerBase) -> int:
765 """
766 Place CopperBottom directly under the given base_layer.
768 base_layer can be:
769 - SubstrateLayer
770 - SilverBottomLayer
771 """
772 if base_layer.surface_tag is None:
773 raise RuntimeError("Base layer must be built before CopperBottom.")
775 if isinstance(base_layer, SubstrateLayer):
776 w = base_layer.width
777 cx = base_layer.center_x
778 cy_base = base_layer.center_y
779 t_base = base_layer.substrate_thickness
781 elif isinstance(base_layer, SilverBottomLayer):
782 w = base_layer.width
783 cx = base_layer.center_x
784 cy_base = base_layer.center_y
785 t_base = base_layer.thickness
787 else:
788 raise TypeError(
789 "CopperBottomLayer.build_under() expected SubstrateLayer or "
790 f"SilverBottomLayer, got {type(base_layer)!r}."
791 )
793 if w is None or cx is None or cy_base is None:
794 raise RuntimeError("Base layer geometry not yet initialized.")
796 # Center of CopperBottom: just below the base layer
797 cy = cy_base - t_base / 2.0 - self.thickness / 2.0
799 self.width = w
800 self.center_x = cx
801 self.center_y = cy
803 x0 = cx - w / 2.0
804 y0 = cy - self.thickness / 2.0
806 self.surface_tag = gmsh.model.occ.addRectangle(x0, y0, 0.0, w, self.thickness)
808 gmsh.model.occ.synchronize()
809 self.refresh_topology()
810 return self.surface_tag
812 def refresh_topology(self) -> None:
813 if self.width is None or self.center_x is None or self.center_y is None:
814 raise RuntimeError("CopperBottom geometry not yet initialized.")
815 super().refresh_topology(
816 self.width,
817 self.thickness,
818 self.center_x,
819 self.center_y,
820 )
822 def create_physical_groups(self, name_prefix: str = "CopperBottom") -> Dict[str, int]:
823 """
824 Create physical groups for the lower copper layer and its edges.
825 """
826 if self.surface_tag is None:
827 raise RuntimeError("build_under() must be called before create_physical_groups().")
829 out: Dict[str, int] = {}
831 pg = gmsh.model.addPhysicalGroup(2, [self.surface_tag])
832 gmsh.model.setPhysicalName(2, pg, name_prefix)
833 out[name_prefix] = pg
835 for name, tags in self.edge_tags.items():
836 edge_pg = gmsh.model.addPhysicalGroup(1, tags)
837 gmsh.model.setPhysicalName(1, edge_pg, f"{name_prefix}_{name}")
838 out[f"{name_prefix}_{name}"] = edge_pg
840 return out
843class CopperTopLayer(_RectLayerBase):
844 """
845 Top 2D copper layer placed directly above:
847 - SilverTop layer if present, otherwise
848 - directly above the HTS layer.
849 """
851 def __init__(self, thickness: float) -> None:
852 super().__init__()
853 self.thickness = float(thickness)
854 self.width: Optional[float] = None
855 self.center_x: Optional[float] = None
856 self.center_y: Optional[float] = None
858 def build_over(self, base_layer: _RectLayerBase) -> int:
859 """
860 Place CopperTop directly above `base_layer`.
862 base_layer can be:
863 - HTSLayer
864 - SilverTopLayer
865 """
866 if base_layer.surface_tag is None:
867 raise RuntimeError("Base layer must be built before CopperTop.")
869 # Case 1: HTS as base
870 if isinstance(base_layer, HTSLayer):
871 w = base_layer.HTS_width
872 cx = base_layer.HTS_center_x
873 cy_base = base_layer.HTS_center_y
874 t_base = base_layer.HTS_thickness
876 # Case 2: SilverTop as base
877 elif isinstance(base_layer, SilverTopLayer):
878 w = base_layer.width
879 cx = base_layer.center_x
880 cy_base = base_layer.center_y
881 t_base = base_layer.thickness
883 else:
884 raise TypeError(
885 "CopperTopLayer.build_over() expected HTSLayer or "
886 f"SilverTopLayer, got {type(base_layer)!r}."
887 )
889 if w is None or cx is None or cy_base is None:
890 raise RuntimeError("Base layer geometry not yet initialized.")
892 cy = cy_base + t_base / 2.0 + self.thickness / 2.0
894 self.width = w
895 self.center_x = cx
896 self.center_y = cy
898 x0 = cx - w / 2.0
899 y0 = cy - self.thickness / 2.0
901 self.surface_tag = gmsh.model.occ.addRectangle(x0, y0, 0.0, w, self.thickness)
903 gmsh.model.occ.synchronize()
904 self.refresh_topology()
905 return self.surface_tag
907 def refresh_topology(self) -> None:
908 if self.width is None or self.center_x is None or self.center_y is None:
909 raise RuntimeError("CopperTop geometry not yet initialized.")
910 super().refresh_topology(
911 self.width,
912 self.thickness,
913 self.center_x,
914 self.center_y,
915 )
917 def create_physical_groups(self, name_prefix: str = "CopperTop") -> Dict[str, int]:
918 """
919 Create physical groups for the top copper layer and its edges.
920 """
921 if self.surface_tag is None:
922 raise RuntimeError("build_over() must be called before create_physical_groups().")
924 out: Dict[str, int] = {}
926 pg = gmsh.model.addPhysicalGroup(2, [self.surface_tag])
927 gmsh.model.setPhysicalName(2, pg, name_prefix)
928 out[name_prefix] = pg
930 for name, tags in self.edge_tags.items():
931 edge_pg = gmsh.model.addPhysicalGroup(1, tags)
932 gmsh.model.setPhysicalName(1, edge_pg, f"{name_prefix}_{name}")
933 out[f"{name_prefix}_{name}"] = edge_pg
935 return out
938class CopperLeftLayer(_RectLayerBase):
939 """
940 Left copper shim, placed directly to the left of the HTS + Substrate
941 (and optionally CopperBottom/CopperTop if present).
943 Height is from:
944 bottom of (CopperBottom or Substrate)
945 to
946 top of (CopperTop or HTS)
947 """
949 def __init__(self, thickness: float) -> None:
950 super().__init__()
951 self.thickness = float(thickness)
952 self.width = self.thickness # Horizontal size
953 self.height: Optional[float] = None
954 self.center_x: Optional[float] = None
955 self.center_y: Optional[float] = None
957 def build_left_of(
958 self,
959 hts: HTSLayer,
960 sub: SubstrateLayer,
961 cu_bottom: Optional["CopperBottomLayer"] = None,
962 cu_top: Optional["CopperTopLayer"] = None,
963 ) -> int:
964 """
965 Build CopperLeft against the left edges of HTS and Substrate,
966 extending vertically to cover the full bottom/top stack.
967 """
968 if hts.surface_tag is None or sub.surface_tag is None:
969 raise RuntimeError("HTS and Substrate must be built before CopperLeft.")
971 if sub.width is None or sub.center_x is None or sub.center_y is None:
972 raise RuntimeError("Substrate geometry not yet initialized.")
974 # Vertical limits
975 if cu_top is not None:
976 y_top = cu_top.center_y + cu_top.thickness / 2.0
977 else:
978 y_top = hts.HTS_center_y + hts.HTS_thickness / 2.0
980 if cu_bottom is not None:
981 y_bot = cu_bottom.center_y - cu_bottom.thickness / 2.0
982 else:
983 y_bot = sub.center_y - sub.substrate_thickness / 2.0
985 self.height = y_top - y_bot
986 if self.height <= 0.0:
987 raise RuntimeError(f"CopperLeft height <= 0 (y_top={y_top}, y_bot={y_bot}).")
989 # Horizontal placement: flush with leftmost face of HTS/Substrate
990 x_left_face = min(
991 hts.HTS_center_x - hts.HTS_width / 2.0,
992 sub.center_x - sub.width / 2.0,
993 )
995 self.center_x = x_left_face - self.thickness / 2.0
996 self.center_y = 0.5 * (y_top + y_bot)
998 x0 = self.center_x - self.thickness / 2.0
999 y0 = self.center_y - self.height / 2.0
1001 self.surface_tag = gmsh.model.occ.addRectangle(x0, y0, 0.0, self.thickness, self.height)
1003 gmsh.model.occ.synchronize()
1004 self._refresh_from_surface()
1005 self._classify_edges_from_bbox(
1006 self.thickness,
1007 self.height,
1008 self.center_x,
1009 self.center_y,
1010 )
1012 return self.surface_tag
1014 def refresh_topology(self) -> None:
1015 if (
1016 self.height is None
1017 or self.center_x is None
1018 or self.center_y is None
1019 ):
1020 raise RuntimeError("CopperLeft geometry not yet initialized.")
1021 self._refresh_from_surface()
1022 self._classify_edges_from_bbox(
1023 self.thickness,
1024 self.height,
1025 self.center_x,
1026 self.center_y,
1027 )
1029 def create_physical_groups(self, name_prefix: str = "CopperLeft") -> Dict[str, int]:
1030 """
1031 Create physical groups for the left copper shim and its edges.
1032 """
1033 if self.surface_tag is None:
1034 raise RuntimeError("build_left_of() must be called before creating PGs.")
1036 out: Dict[str, int] = {}
1038 pg = gmsh.model.addPhysicalGroup(2, [self.surface_tag])
1039 gmsh.model.setPhysicalName(2, pg, name_prefix)
1040 out[name_prefix] = pg
1042 for name, tags in self.edge_tags.items():
1043 edge_pg = gmsh.model.addPhysicalGroup(1, tags)
1044 gmsh.model.setPhysicalName(1, edge_pg, f"{name_prefix}_{name}")
1045 out[f"{name_prefix}_{name}"] = edge_pg
1047 return out
1050class CopperRightLayer(_RectLayerBase):
1051 """
1052 Right copper shim, symmetric counterpart of CopperLeftLayer.
1053 """
1055 def __init__(self, thickness: float) -> None:
1056 super().__init__()
1057 self.thickness = float(thickness)
1058 self.width = self.thickness
1059 self.height: Optional[float] = None
1060 self.center_x: Optional[float] = None
1061 self.center_y: Optional[float] = None
1063 def build_right_of(
1064 self,
1065 hts: HTSLayer,
1066 sub: SubstrateLayer,
1067 cu_bottom: Optional["CopperBottomLayer"] = None,
1068 cu_top: Optional["CopperTopLayer"] = None,
1069 ) -> int:
1070 """
1071 Build CopperRight against the right edges of HTS and Substrate.
1072 """
1073 if hts.surface_tag is None or sub.surface_tag is None:
1074 raise RuntimeError("HTS and Substrate must be built before CopperRight.")
1076 if sub.width is None or sub.center_x is None or sub.center_y is None:
1077 raise RuntimeError("Substrate geometry not yet initialized.")
1079 if cu_top is not None:
1080 y_top = cu_top.center_y + cu_top.thickness / 2.0
1081 else:
1082 y_top = hts.HTS_center_y + hts.HTS_thickness / 2.0
1084 if cu_bottom is not None:
1085 y_bot = cu_bottom.center_y - cu_bottom.thickness / 2.0
1086 else:
1087 y_bot = sub.center_y - sub.substrate_thickness / 2.0
1089 self.height = y_top - y_bot
1090 if self.height <= 0.0:
1091 raise RuntimeError(f"CopperRight height <= 0 (y_top={y_top}, y_bot={y_bot}).")
1093 x_right_face = max(
1094 hts.HTS_center_x + hts.HTS_width / 2.0,
1095 sub.center_x + sub.width / 2.0,
1096 )
1098 self.center_x = x_right_face + self.thickness / 2.0
1099 self.center_y = 0.5 * (y_top + y_bot)
1101 x0 = self.center_x - self.thickness / 2.0
1102 y0 = self.center_y - self.height / 2.0
1104 self.surface_tag = gmsh.model.occ.addRectangle(x0, y0, 0.0, self.thickness, self.height)
1106 gmsh.model.occ.synchronize()
1107 self._refresh_from_surface()
1108 self._classify_edges_from_bbox(
1109 self.thickness,
1110 self.height,
1111 self.center_x,
1112 self.center_y,
1113 )
1115 return self.surface_tag
1117 def refresh_topology(self) -> None:
1118 if (
1119 self.height is None
1120 or self.center_x is None
1121 or self.center_y is None
1122 ):
1123 raise RuntimeError("CopperRight geometry not yet initialized.")
1124 self._refresh_from_surface()
1125 self._classify_edges_from_bbox(
1126 self.thickness,
1127 self.height,
1128 self.center_x,
1129 self.center_y,
1130 )
1132 def create_physical_groups(self, name_prefix: str = "CopperRight") -> Dict[str, int]:
1133 """
1134 Create physical groups for the right copper shim and its edges.
1135 """
1136 if self.surface_tag is None:
1137 raise RuntimeError("build_right_of() must be called before creating PGs.")
1139 out: Dict[str, int] = {}
1141 pg = gmsh.model.addPhysicalGroup(2, [self.surface_tag])
1142 gmsh.model.setPhysicalName(2, pg, name_prefix)
1143 out[name_prefix] = pg
1145 for name, tags in self.edge_tags.items():
1146 edge_pg = gmsh.model.addPhysicalGroup(1, tags)
1147 gmsh.model.setPhysicalName(1, edge_pg, f"{name_prefix}_{name}")
1148 out[f"{name_prefix}_{name}"] = edge_pg
1150 return out
1153# ---------------------------------------------------------------------------
1154# High-level geometry builder
1155# ---------------------------------------------------------------------------
1157class Generate_geometry:
1158 """
1159 Generate a 2D geometry of the CACCC model and save as .brep / .xao.
1161 Responsibilities
1162 ----------------
1163 - Initialize a Gmsh OCC model for the given magnet_name.
1164 - Instantiate and build each layer:
1165 HTS, Substrate, optional SilverTop / SilverBottom,
1166 optional CopperTop / CopperBottom / CopperLeft / CopperRight.
1167 - Create an enclosing circular air region and cut out all solids.
1168 - Refresh topology and create physical groups:
1169 - material regions (2D)
1170 - layer edges (1D)
1171 - inter-layer interfaces (1D)
1172 - air outer boundary (1D)
1173 - Write:
1174 <magnet_name>.brep
1175 <magnet_name>.xao
1176 """
1178 def __init__(
1179 self,
1180 fdm: FDM,
1181 inputs_folder_path: str,
1182 verbose: bool = True,
1183 *,
1184 initialize_gmsh: bool = True,
1185 create_model: bool = True,
1186 create_physical_groups: Optional[bool] = None,
1187 wipe_physical_groups: Optional[bool] = None,
1188 write_files: Optional[bool] = None,
1189 clear_gmsh_on_finalize: Optional[bool] = None,
1190 external_gu: Optional[GmshUtils] = None,
1191 ) -> None:
1193 self.fdm = fdm
1194 self.inputs_folder_path = inputs_folder_path
1196 self.model_folder = os.path.join(os.getcwd())
1197 self.magnet_name = fdm.general.magnet_name
1199 self.model_file = os.path.join(self.model_folder, f"{self.magnet_name}.brep")
1200 self.xao_file = os.path.join(self.model_folder, f"{self.magnet_name}.xao")
1202 self.verbose = verbose
1203 self.gu = external_gu if external_gu is not None else GmshUtils(self.model_folder, self.verbose)
1205 # When embedded in another builder (e.g. TSTC 2D stack),
1206 # we must not re-initialize Gmsh, and must not create a new model.
1207 self._create_model = bool(create_model)
1210 embedded = not self._create_model
1212 # Standalone defaults (embedded=False): behave like before.
1213 # Embedded defaults (embedded=True): do NOT clear / wipe / write unless requested.
1214 self._create_physical_groups = (
1215 bool(create_physical_groups)
1216 if create_physical_groups is not None
1217 else (not embedded)
1218 )
1219 self._wipe_physical_groups = (
1220 bool(wipe_physical_groups)
1221 if wipe_physical_groups is not None
1222 else (not embedded)
1223 )
1224 self._write_files = (
1225 bool(write_files)
1226 if write_files is not None
1227 else (not embedded)
1228 )
1229 self._clear_gmsh_on_finalize = (
1230 bool(clear_gmsh_on_finalize)
1231 if clear_gmsh_on_finalize is not None
1232 else (not embedded)
1233 )
1236 conductor_dict = self.fdm.conductors
1238 if not conductor_dict:
1239 raise KeyError(
1240 "fdm.conductors is empty: cannot build CC geometry. "
1241 "Expected at least one conductor entry."
1242 )
1244 # Preferred legacy key (many FiQuS models use this)
1245 selected_conductor_name = self.fdm.magnet.solve.conductor_name
1246 if selected_conductor_name in conductor_dict:
1247 selected_conductor: Conductor = conductor_dict[selected_conductor_name]
1248 else:
1249 raise ValueError(
1250 f"Conductor name: {selected_conductor_name} not present in the conductors section"
1251 )
1253 if initialize_gmsh:
1254 self.gu.initialize(verbosity_Gmsh=fdm.run.verbosity_Gmsh)
1256 # OCC warnings (0 = silenced in terminal)
1257 gmsh.option.setNumber("General.Terminal", 0)
1258 gmsh.option.setNumber("General.Verbosity", 0)
1260 # All built layers are stored here (by name)
1261 self.layers: Dict[str, object] = {}
1262 self._model_ready = False
1265 # Reference center for this CC cross-section (used for placement and air centering)
1266 self._geom_center_x: float = 0.0
1267 self._geom_center_y: float = 0.0
1270 # the strand is a Union[Round, Rectangular, CC, Homogenized]
1271 strand = selected_conductor.strand
1272 if not isinstance(strand, CC):
1273 raise TypeError(
1274 f"Expected strand type 'CC' for CACCC geometry, got {type(strand)}"
1275 )
1277 # Store this CC strand so all generate_* methods can use it
1278 self.cc_strand = strand
1279 s = self.cc_strand # shorthand
1281 # ------------------------------------------------------------------
1282 # Geometry sanity checks (enforced here, not in magnet.geometry model)
1283 # ------------------------------------------------------------------
1284 # Mandatory base layer dimensions must be > 0
1285 _require_positive(
1286 s.HTS_width,
1287 name=f"conductors.{selected_conductor_name}.strand.HTS_width",
1288 )
1289 _require_positive(
1290 s.HTS_thickness,
1291 name=f"conductors.{selected_conductor_name}.strand.HTS_thickness",
1292 )
1293 _require_positive(
1294 s.substrate_thickness,
1295 name=f"conductors.{selected_conductor_name}.strand.substrate_thickness",
1296 )
1298 # Striation parameters
1299 n_fil = _require_int_ge(
1300 getattr(s, "number_of_filaments", None),
1301 name=f"conductors.{selected_conductor_name}.strand.number_of_filaments",
1302 min_value=1,
1303 )
1304 gap = getattr(s, "gap_between_filaments", None)
1305 gap_v = _require_non_negative_optional(
1306 gap,
1307 name=f"conductors.{selected_conductor_name}.strand.gap_between_filaments",
1308 )
1310 if n_fil is not None and n_fil > 1:
1311 if gap is None or gap_v <= 0.0:
1312 raise ValueError(
1313 "HTS striation requested (number_of_filaments > 1) but "
1314 "gap_between_filaments is missing or <= 0."
1315 )
1317 total_groove = gap_v * (n_fil - 1)
1318 hts_w = float(s.HTS_width)
1319 if total_groove >= hts_w:
1320 raise ValueError(
1321 f"Total groove width {total_groove:g} >= HTS_width {hts_w:g}. "
1322 "Reduce gap_between_filaments or number_of_filaments."
1323 )
1325 # ------------------------------------------------------------------
1326 # Layer inclusion logic (thickness-based only, but negatives are errors)
1327 # ------------------------------------------------------------------
1328 ag_top = _require_non_negative_optional(
1329 getattr(getattr(s, "silver_thickness", None), "top", None),
1330 name=f"conductors.{selected_conductor_name}.strand.silver_thickness.top",
1331 )
1332 ag_bot = _require_non_negative_optional(
1333 getattr(getattr(s, "silver_thickness", None), "bottom", None),
1334 name=f"conductors.{selected_conductor_name}.strand.silver_thickness.bottom",
1335 )
1337 cu_top = _require_non_negative_optional(
1338 getattr(getattr(s, "copper_thickness", None), "top", None),
1339 name=f"conductors.{selected_conductor_name}.strand.copper_thickness.top",
1340 )
1341 cu_bot = _require_non_negative_optional(
1342 getattr(getattr(s, "copper_thickness", None), "bottom", None),
1343 name=f"conductors.{selected_conductor_name}.strand.copper_thickness.bottom",
1344 )
1345 cu_left = _require_non_negative_optional(
1346 getattr(getattr(s, "copper_thickness", None), "left", None),
1347 name=f"conductors.{selected_conductor_name}.strand.copper_thickness.left",
1348 )
1349 cu_right = _require_non_negative_optional(
1350 getattr(getattr(s, "copper_thickness", None), "right", None),
1351 name=f"conductors.{selected_conductor_name}.strand.copper_thickness.right",
1352 )
1354 self._use_silver_top = ag_top > 0.0
1355 self._use_silver_bottom = ag_bot > 0.0
1357 self._use_copper_top = cu_top > 0.0
1358 self._use_copper_bottom = cu_bot > 0.0
1359 self._use_copper_left = cu_left > 0.0
1360 self._use_copper_right = cu_right > 0.0
1365 # Internal helpers
1366 def _ensure_model(self) -> None:
1367 """
1368 Ensure there is an active Gmsh model.
1370 In embedded mode (create_model=False), we assume the caller already did:
1371 gmsh.model.add(...)
1372 """
1373 if self._model_ready:
1374 return
1376 if getattr(self, "_create_model", True):
1377 gmsh.model.add(self.magnet_name)
1379 self._model_ready = True
1382 # Layer builders
1383 def generate_HTS_layer(self, center_x: float = 0.0, center_y: float = 0.0) -> None:
1384 """
1385 Build a single HTS layer (2D) centered at (center_x, center_y).
1387 Striation parameters (from CC strand):
1388 - number_of_filaments
1389 - gap_between_filaments
1390 """
1391 self._ensure_model()
1393 # Save reference center so other steps (e.g. air disk) can follow placement.
1394 self._geom_center_x = float(center_x)
1395 self._geom_center_y = float(center_y)
1397 s = self.cc_strand # CC instance
1399 layer = HTSLayer(
1400 HTS_thickness=float(s.HTS_thickness),
1401 HTS_width=float(s.HTS_width),
1402 HTS_center_x=float(center_x),
1403 HTS_center_y=float(center_y),
1404 number_of_filaments=s.number_of_filaments,
1405 gap_between_filaments=s.gap_between_filaments,
1406 )
1408 layer.build_HTS()
1409 self.layers["HTS"] = layer
1414 def generate_silver_top_layer(self) -> None:
1415 """
1416 Build the SilverTop layer directly above HTS (if enabled).
1417 """
1418 self._ensure_model()
1420 if not self._use_silver_top:
1421 logger.info("[Geom-CC] Skipping SilverTop layer (thickness == 0).")
1422 return
1424 if "HTS" not in self.layers:
1425 raise RuntimeError("HTS must be built before top Silver.")
1427 s = self.cc_strand
1429 Ag = SilverTopLayer(thickness=float(s.silver_thickness.top))
1430 Ag.build_over(self.layers["HTS"])
1431 self.layers["SilverTop"] = Ag
1435 def generate_substrate_layer(self) -> None:
1436 """
1437 Build the substrate layer directly under the HTS (shared interface).
1438 """
1439 self._ensure_model()
1441 if "HTS" not in self.layers:
1442 raise RuntimeError(
1443 "HTS layer must be built before substrate. "
1444 "Call generate_HTS_layer() first."
1445 )
1447 hts = self.layers["HTS"]
1448 s = self.cc_strand
1450 substrate = SubstrateLayer(substrate_thickness=float(s.substrate_thickness))
1451 substrate.build_substrate(hts)
1452 self.layers["Substrate"] = substrate
1455 def generate_silver_bottom_layer(self) -> None:
1456 """
1457 Build the bottom Silver layer directly under the Substrate
1458 (shared interface), if enabled.
1459 """
1460 self._ensure_model()
1462 if not self._use_silver_bottom:
1463 logger.info("[Geom-CC] Skipping SilverBottom layer (config disabled).")
1464 return
1466 if "Substrate" not in self.layers:
1467 raise RuntimeError("Substrate must be built before bottom Silver.")
1469 sub = self.layers["Substrate"]
1470 s = self.cc_strand
1472 AgB = SilverBottomLayer(thickness=float(s.silver_thickness.bottom))
1473 AgB.build_under(sub)
1474 self.layers["SilverBottom"] = AgB
1478 def generate_copper_bottom_layer(self) -> None:
1479 """
1480 Build CopperBottom directly under the bottom-most central layer:
1482 - Under SilverBottom if it exists,
1483 - Otherwise directly under Substrate.
1485 (Side shims are handled separately.)
1486 """
1487 self._ensure_model()
1489 if not self._use_copper_bottom:
1490 logger.info("[Geom-CC] Skipping CopperBottom layer (config disabled).")
1491 return
1493 s = self.cc_strand
1495 if "Substrate" not in self.layers:
1496 raise RuntimeError("Substrate must be built before CopperBottom.")
1498 base_layer = self.layers.get("SilverBottom", None)
1499 if base_layer is None:
1500 base_layer = self.layers["Substrate"]
1502 cuL = CopperBottomLayer(thickness=float(s.copper_thickness.bottom))
1503 cuL.build_under(base_layer)
1504 self.layers["CopperBottom"] = cuL
1507 def generate_copper_top_layer(self) -> None:
1508 """
1509 Build CopperTop directly above HTS (or above SilverTop if present).
1510 """
1511 self._ensure_model()
1513 if not self._use_copper_top:
1514 logger.info("[Geom-CC] Skipping CopperTop layer (config disabled).")
1515 return
1517 s = self.cc_strand
1519 if "HTS" not in self.layers:
1520 raise RuntimeError("HTS must be built before CopperTop.")
1522 hts = self.layers["HTS"]
1523 cuT = CopperTopLayer(thickness=float(s.copper_thickness.top))
1524 silver_top = self.layers.get("SilverTop", None)
1526 if silver_top is not None:
1527 cuT.build_over(silver_top)
1528 else:
1529 cuT.build_over(hts)
1531 self.layers["CopperTop"] = cuT
1535 def generate_copper_left_layer(self) -> None:
1536 """
1537 Build CopperLeft against the left edges of HTS + Substrate,
1538 extended vertically to cover the whole stack (CuBottom / CuTop if present).
1539 """
1540 self._ensure_model()
1542 if not self._use_copper_left:
1543 logger.info("[Geom-CC] Skipping CopperLeft layer (config disabled).")
1544 return
1546 s = self.cc_strand
1548 if "HTS" not in self.layers or "Substrate" not in self.layers:
1549 raise RuntimeError("Build HTS and Substrate before CopperLeft.")
1551 cuL = CopperLeftLayer(thickness=float(s.copper_thickness.left))
1553 cuL.build_left_of(
1554 self.layers["HTS"],
1555 self.layers["Substrate"],
1556 cu_bottom=self.layers.get("CopperBottom"),
1557 cu_top=self.layers.get("CopperTop"),
1558 )
1560 self.layers["CopperLeft"] = cuL
1563 def generate_copper_right_layer(self) -> None:
1564 """
1565 Build CopperRight against the right edges of HTS + Substrate,
1566 extended vertically to cover the whole stack (CuBottom / CuTop if present).
1567 """
1568 self._ensure_model()
1570 if not self._use_copper_right:
1571 logger.info("[Geom-CC] Skipping CopperRight layer (config disabled).")
1572 return
1574 s = self.cc_strand
1576 if "HTS" not in self.layers or "Substrate" not in self.layers:
1577 raise RuntimeError("Build HTS and Substrate before CopperRight.")
1579 cuR = CopperRightLayer(thickness=float(s.copper_thickness.right))
1581 cuR.build_right_of(
1582 self.layers["HTS"],
1583 self.layers["Substrate"],
1584 cu_bottom=self.layers.get("CopperBottom"),
1585 cu_top=self.layers.get("CopperTop"),
1586 )
1588 self.layers["CopperRight"] = cuR
1591 # Air region
1593 def generate_air_region(self) -> None:
1594 g = self.fdm.magnet.geometry
1595 R = _require_positive(g.air_radius, name="magnet.geometry.air_radius")
1597 # -------------------------------
1598 # 1) Collect all inner surfaces
1599 # -------------------------------
1600 inner: List[int] = []
1601 for L in self.layers.values():
1602 if hasattr(L, "filament_tags") and getattr(L, "filament_tags"):
1603 inner.extend(getattr(L, "filament_tags"))
1605 if hasattr(L, "groove_tags") and getattr(L, "groove_tags"):
1606 inner.extend(getattr(L, "groove_tags"))
1608 elif getattr(L, "surface_tag", None) is not None:
1609 tag = getattr(L, "surface_tag")
1611 if tag not in inner:
1612 inner.append(tag)
1614 # Defensive: collapse exact duplicates before any boolean ops
1615 gmsh.model.occ.removeAllDuplicates()
1616 gmsh.model.occ.synchronize()
1618 # -----------------------------------------------------
1619 # 2) Enforce conformity inside the conductor stack
1620 # (split all touching rectangles so they share edges)
1621 # -----------------------------------------------------
1622 # Fragmenting the conductor stack first enforces conformal interfaces:
1623 # touching layers/filaments share identical curve entities, which makes
1624 # interface physical groups and meshing constraints reliable.
1625 all_rects = [(2, s) for s in inner]
1626 frag_conduct, _ = gmsh.model.occ.fragment(all_rects, [])
1627 gmsh.model.occ.synchronize()
1630 # Rebuild 'inner' from the fragment result (their tags can change)
1631 inner = [ent[1] for ent in frag_conduct]
1633 # ------------------------------------
1634 # 3) Create the Air disk
1635 # ------------------------------------
1636 # Robust air centering:
1637 # compute the bbox of the *actual* conductor surfaces and center the air on that.
1638 xmin = float("inf")
1639 ymin = float("inf")
1640 xmax = float("-inf")
1641 ymax = float("-inf")
1643 for s_tag in inner:
1644 bxmin, bymin, _, bxmax, bymax, _ = gmsh.model.getBoundingBox(2, int(s_tag))
1645 xmin = min(xmin, float(bxmin))
1646 ymin = min(ymin, float(bymin))
1647 xmax = max(xmax, float(bxmax))
1648 ymax = max(ymax, float(bymax))
1650 if not (xmin < xmax and ymin < ymax):
1651 # Fallback (should not happen unless inner is empty / invalid)
1652 cx = float(getattr(self, "_geom_center_x", 0.0))
1653 cy = float(getattr(self, "_geom_center_y", 0.0))
1654 else:
1655 cx = 0.5 * (xmin + xmax)
1656 cy = 0.5 * (ymin + ymax)
1658 disk = gmsh.model.occ.addDisk(cx, cy, 0.0, R, R)
1661 # Fragment(disk, inner) splits the disk into multiple surfaces:
1662 # one is the exterior air region, others may be trapped pockets.
1663 # The largest disk-derived surface is taken as the physical Air domain.
1664 out_dimtags, out_map = gmsh.model.occ.fragment(
1665 [(2, disk)],
1666 [(2, s) for s in inner],
1667 )
1668 gmsh.model.occ.synchronize()
1671 # out_map[0] corresponds to the object list = [(2, disk)]
1672 disk_pieces = out_map[0] if out_map and out_map[0] else []
1673 if not disk_pieces:
1674 raise RuntimeError("Air fragment produced no disk-derived pieces (unexpected).")
1676 # Pick the largest disk-derived piece as the Air region
1677 air_ent = max(
1678 disk_pieces,
1679 key=lambda e: float(gmsh.model.occ.getMass(e[0], e[1])),
1680 )
1681 air_tag = int(air_ent[1])
1683 self.layers["Air"] = SimpleNamespace(surface_tag=air_tag)
1687 # Topology + physical groups
1689 def _refresh_all(self) -> None:
1690 """
1691 After OCC cleanup, update curve/point lists and edge classification
1692 for all rectangular layers that implement refresh_topology().
1693 """
1694 for L in self.layers.values():
1695 if hasattr(L, "refresh_topology"):
1696 L.refresh_topology() # type: ignore[call-arg]
1699 def _create_all_physical_groups(self, name_prefix: str = "") -> None:
1700 """
1701 Create all material + interface + air physical groups.
1703 If name_prefix != "":
1704 all physical group names become: "<name_prefix>_<OriginalName>"
1705 and spaces in the original name are replaced by underscores.
1706 """
1707 prefix = str(name_prefix).strip()
1709 # Embedded mode: do not wipe physical groups by default, because other builders
1710 # may have already created groups in the same Gmsh model.
1711 # Use name_prefix to avoid naming collisions across repeated CC instances.
1712 if self._wipe_physical_groups:
1713 gmsh.model.removePhysicalGroups()
1714 elif not prefix:
1715 raise ValueError(
1716 "name_prefix must be non-empty when wipe_physical_groups=False "
1717 "(otherwise physical group names will collide)."
1718 )
1722 prefix = str(name_prefix).strip()
1724 def _pref(name: str) -> str:
1725 if not prefix:
1726 return name
1727 safe = str(name).replace(" ", "_")
1728 return f"{prefix}_{safe}"
1730 # Layer physical groups (2D + edges)
1731 if "HTS" in self.layers:
1732 self.layers["HTS"].create_physical_groups(_pref("HTS")) # type: ignore[arg-type]
1734 if "Substrate" in self.layers:
1735 self.layers["Substrate"].create_physical_groups(_pref("Substrate")) # type: ignore[arg-type]
1737 if "SilverTop" in self.layers:
1738 self.layers["SilverTop"].create_physical_groups(_pref("SilverTop")) # type: ignore[arg-type]
1740 if "SilverBottom" in self.layers:
1741 self.layers["SilverBottom"].create_physical_groups(_pref("SilverBottom")) # type: ignore[arg-type]
1743 if "CopperBottom" in self.layers:
1744 self.layers["CopperBottom"].create_physical_groups(_pref("CopperBottom")) # type: ignore[arg-type]
1746 if "CopperTop" in self.layers:
1747 self.layers["CopperTop"].create_physical_groups(_pref("CopperTop")) # type: ignore[arg-type]
1749 if "CopperLeft" in self.layers:
1750 self.layers["CopperLeft"].create_physical_groups(_pref("CopperLeft")) # type: ignore[arg-type]
1752 if "CopperRight" in self.layers:
1753 self.layers["CopperRight"].create_physical_groups(_pref("CopperRight")) # type: ignore[arg-type]
1755 # Interfaces (1D)
1756 # Interface convention:
1757 # We define each interface using the boundary curves of ONE of the adjacent
1758 # layers (typically the "lower" layer's upper edge or vice versa). Because we
1759 # fragment the stack earlier, both sides should reference the same curve IDs.
1760 # Interfaces (1D)
1762 if "Substrate" in self.layers and "SilverBottom" in self.layers:
1763 iface = self.layers["Substrate"].edge_tags["Lower"] # type: ignore[index]
1764 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1765 gmsh.model.setPhysicalName(1, pg_if, _pref("Substrate_SilverBottom_Interface"))
1767 if "SilverBottom" in self.layers and "CopperBottom" in self.layers:
1768 iface = self.layers["SilverBottom"].edge_tags["Lower"] # type: ignore[index]
1769 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1770 gmsh.model.setPhysicalName(1, pg_if, _pref("SilverBottom_CopperBottom_Interface"))
1772 if (
1773 "Substrate" in self.layers
1774 and "CopperBottom" in self.layers
1775 and "SilverBottom" not in self.layers
1776 ):
1777 iface = self.layers["CopperBottom"].edge_tags["Upper"] # type: ignore[index]
1778 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1779 gmsh.model.setPhysicalName(1, pg_if, _pref("Substrate_CopperBottom_Interface"))
1781 if "HTS" in self.layers and "Substrate" in self.layers:
1782 iface = self.layers["Substrate"].edge_tags["Upper"] # type: ignore[index]
1783 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1784 gmsh.model.setPhysicalName(1, pg_if, _pref("Buffer Layer"))
1786 if "HTS" in self.layers and "SilverTop" in self.layers:
1787 iface = self.layers["HTS"].edge_tags["Upper"] # type: ignore[index]
1788 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1789 gmsh.model.setPhysicalName(1, pg_if, _pref("HTS_Silver_Interface"))
1791 if "SilverTop" in self.layers and "CopperTop" in self.layers:
1792 iface = self.layers["SilverTop"].edge_tags["Upper"] # type: ignore[index]
1793 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1794 gmsh.model.setPhysicalName(1, pg_if, _pref("SilverTop_CopperTop_Interface"))
1796 if (
1797 "HTS" in self.layers
1798 and "CopperTop" in self.layers
1799 and "SilverTop" not in self.layers
1800 ):
1801 iface = self.layers["HTS"].edge_tags["Upper"] # type: ignore[index]
1802 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1803 gmsh.model.setPhysicalName(1, pg_if, _pref("HTS_CopperTop_Interface"))
1805 if "CopperLeft" in self.layers and "HTS" in self.layers:
1806 iface = self.layers["HTS"].edge_tags["LeftEdge"] # type: ignore[index]
1807 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1808 gmsh.model.setPhysicalName(1, pg_if, _pref("HTS_CopperLeft_Interface"))
1810 if "CopperLeft" in self.layers and "Substrate" in self.layers:
1811 iface = self.layers["Substrate"].edge_tags["LeftEdge"] # type: ignore[index]
1812 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1813 gmsh.model.setPhysicalName(1, pg_if, _pref("Substrate_CopperLeft_Interface"))
1815 if "CopperRight" in self.layers and "HTS" in self.layers:
1816 iface = self.layers["HTS"].edge_tags["RightEdge"] # type: ignore[index]
1817 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1818 gmsh.model.setPhysicalName(1, pg_if, _pref("HTS_CopperRight_Interface"))
1820 if "CopperRight" in self.layers and "Substrate" in self.layers:
1821 iface = self.layers["Substrate"].edge_tags["RightEdge"] # type: ignore[index]
1822 pg_if = gmsh.model.addPhysicalGroup(1, iface)
1823 gmsh.model.setPhysicalName(1, pg_if, _pref("Substrate_CopperRight_Interface"))
1825 # Air (surface + outer/inner boundary + inner boundary + one point for gauging)
1826 if "Air" in self.layers:
1827 air_tag = self.layers["Air"].surface_tag # type: ignore[attr-defined]
1829 pg_s = gmsh.model.addPhysicalGroup(2, [air_tag])
1830 gmsh.model.setPhysicalName(2, pg_s, _pref("Air"))
1832 b = gmsh.model.getBoundary([(2, air_tag)], oriented=False, recursive=False)
1833 cand = [c[1] for c in b]
1835 outer_curves: List[int] = []
1836 for ctag in cand:
1837 try:
1838 _, up = gmsh.model.getAdjacencies(1, ctag)
1839 if len(up) == 1 and up[0] == air_tag:
1840 outer_curves.append(ctag)
1841 except Exception:
1842 pass
1844 if not outer_curves:
1845 lengths = [(ctag, gmsh.model.occ.getMass(1, ctag)) for ctag in cand]
1846 max_len = max((L for _, L in lengths), default=0.0)
1847 if max_len > 0.0:
1848 outer_curves = [
1849 ctag for ctag, L in lengths
1850 if abs(L - max_len) <= 1e-9 * max_len
1851 ]
1853 pg_c = gmsh.model.addPhysicalGroup(1, outer_curves)
1854 gmsh.model.setPhysicalName(1, pg_c, _pref("Air_Outer"))
1856 inner_curves = [ctag for ctag in cand if ctag not in outer_curves]
1857 pg_c = gmsh.model.addPhysicalGroup(1, inner_curves)
1858 gmsh.model.setPhysicalName(1, pg_c, _pref("Air_Inner"))
1861 # One point on the air boundary for gauging the scalar magnetic potential phi
1862 adj = gmsh.model.getAdjacencies(1, outer_curves[0])
1863 point = adj[1][0] # downward, [1], to get dimension 0 and first element to grab a single point tag
1864 pg_c = gmsh.model.addPhysicalGroup(0, [point])
1865 gmsh.model.setPhysicalName(0, pg_c, _pref("Gauging_point"))
1869 # Finalize
1870 def finalize_and_write(self, name_prefix: str = "") -> None:
1871 """
1872 Perform final OCC synchronize, refresh topology, recreate all
1873 physical groups, and write .brep and .xao.
1874 """
1875 gmsh.model.occ.removeAllDuplicates()
1876 gmsh.model.occ.synchronize()
1878 # Refresh topology after final OCC operations.
1879 self._refresh_all()
1881 # Create physical groups BEFORE writing .xao so they are embedded.
1882 if self._create_physical_groups:
1883 self._create_all_physical_groups(name_prefix=name_prefix)
1885 # Only write files if requested (standalone default: True).
1886 if self._write_files:
1887 gmsh.write(self.model_file)
1888 gmsh.write(self.xao_file)
1890 if self.fdm.run.launch_gui and self._create_model:
1891 self.gu.launch_interactive_GUI()
1892 elif self._clear_gmsh_on_finalize:
1893 gmsh.clear()
1899 # Loading an existing geometry
1900 def load_geometry(self, gui: bool = False) -> None:
1901 """
1902 Load an existing .brep geometry and optionally launch the GUI.
1903 """
1904 gmsh.open(self.model_file)
1905 if gui:
1906 self.gu.launch_interactive_GUI()