Coverage for fiqus/mesh_generators/MeshConductorAC_CC.py: 80%
485 statements
« prev ^ index » next coverage.py v7.4.4, created at 2026-01-09 01:49 +0000
1"""
2MeshConductorAC_CC.py
42D mesh generator for the coated-conductor (CACCC) model.
6Assumptions / contract with GeometryConductorAC_CC
7--------------------------------------------------
8- Geometry is axis-aligned: x = tape width, y = tape thickness.
9- The imported .xao file must contain the physical groups expected by this
10 module (e.g. HTS, optional HTS_Grooves, Substrate, SilverTop/Bottom,
11 CopperTop/Bottom/Left/Right, Air, Air_Outer, Air_Inner, and boundary/interface
12 curve groups).
14Meshing strategy
15----------------
16- HTS: structured (transfinite) meshing. If HTS is striated, filaments are
17 meshed as transfinite quads, while grooves can remain unstructured triangles
18 (but with aligned boundary node counts).
19- Copper: uniform target element size via a background field (single knob:
20 copper_elem_scale). We intentionally avoid transfinite on copper to preserve
21 uniform spacing.
22- Air: optionally controlled by air_boundary_mesh_size_ratio on Air_Outer.
23 Boundary sizing can act as a fallback when no background field is used.
24"""
27import os
28import logging
29from typing import Any, Dict, List, Optional
31import gmsh
33from fiqus.data import RegionsModelFiQuS # noqa: F401 (used later for regions)
34from fiqus.utils.Utils import GmshUtils, FilesAndFolders # noqa: F401
35from fiqus.data.RegionsModelFiQuS import RegionsModel # noqa: F401
36from fiqus.data.DataConductor import CC, Conductor
39from fiqus.mesh_generators.MeshConductorAC_Strand import Mesh as BaseMesh
41logger = logging.getLogger("FiQuS")
43occ = gmsh.model.occ
46class Mesh(BaseMesh):
47 """
48 Mesh generator for the CACCC (coated-conductor) 2D cross-section.
49 This class reuses the conductor-AC base Mesh from
50 MeshConductorAC_Strand.
51 """
53 def __init__(self, fdm, verbose: bool = True) -> None:
54 """
55 Initialize the CACCC mesh generator.
57 Parameters
58 ----------
59 fdm : object
60 FiQuS data model instance.
61 verbose : bool, optional
62 If True, gmsh terminal output is enabled. Defaults to True.
63 """
64 super().__init__(fdm, verbose)
66 # Shortcuts to CACCC-specific geometry and mesh sections.
67 # These correspond to DataFiQuSConductorAC_CC.CACCCGeometry
68 # and DataFiQuSConductorAC_CC.CACCCMesh.
69 self.cc_geom = self.cacdm.geometry
70 self.cc_mesh = self.cacdm.mesh
72 # Pull the CC strand from conductors
73 conductor_dict = self.fdm.conductors
74 selected_conductor_name = self.fdm.magnet.solve.conductor_name
75 if selected_conductor_name in conductor_dict:
76 selected_conductor: Conductor = conductor_dict[selected_conductor_name]
77 else:
78 raise ValueError(
79 f"Conductor name: {selected_conductor_name} not present in the conductors section"
80 )
81 strand = selected_conductor.strand
83 if not isinstance(strand, CC):
84 raise TypeError(
85 f"Expected strand type 'CC' for CACCC mesh, got {type(strand)}"
86 )
88 self.cc_strand = strand # CC instance (HTS_width, HTS_thickness, etc.)
90 # Helper dictionaries populated in generate_mesh().
91 # They will be reused in later steps (mesh size fields, regions, ...).
92 self.physical_surfaces_by_name: Dict[str, Dict[str, Any]] = {}
93 self.physical_curves_by_name: Dict[str, Dict[str, Any]] = {}
94 self.physical_points_by_name: Dict[str, Dict[str, Any]] = {}
96 # HTS base mesh sizes (x and y) computed from geometry + mesh parameters.
97 self.hts_base_size_x: Optional[float] = None
98 self.hts_base_size_y: Optional[float] = None
100 # Canonical number of nodes along the tape width (x-direction),
101 # used for HTS and later propagated to other layers.
102 self.horizontal_nodes_x: Optional[int] = None
104 # Equivalent superconducting width (sum of filament widths only)
105 # and corresponding reference horizontal element size.
106 self.equiv_hts_filament_width: Optional[float] = None
107 self.hts_ref_size_x_filaments: Optional[float] = None
109 # Horizontal element densities (elements per unit length) for
110 # HTS filaments and grooves.
111 self.hts_density_x_filaments: Optional[float] = None
112 self.hts_density_x_grooves: Optional[float] = None
114 # Per-surface horizontal node counts (top/bottom HTS edges).
115 # Keys are surface tags (2D entities).
116 self.hts_nodes_x_filaments: Dict[int, int] = {}
117 self.hts_nodes_x_grooves: Dict[int, int] = {}
119 # Approximate vertical discretisation per stacked layer
120 # (Substrate, Silver*, CopperTop/Bottom, HTS). These are
121 # derived from the HTS vertical reference size and the
122 # per-layer *_elem_scale parameters.
123 self.vertical_elems_by_layer: Dict[str, int] = {}
124 self.layer_base_size_y: Dict[str, float] = {}
127 # ------------------------------------------------------------------ #
128 # Internal helpers: physical groups & mesh parameters
129 # ------------------------------------------------------------------ #
131 def _collect_physical_groups(self) -> None:
132 """
133 Populate helper dictionaries for physical groups (surfaces/curves/points).
135 This method scans all physical groups present in the current gmsh model
136 and fills:
137 - self.physical_surfaces_by_name (dim=2)
138 - self.physical_curves_by_name (dim=1)
139 - self.physical_points_by_name (dim=0)
141 Keys are the physical names (strings). Values store both:
142 - phys_tag : used later in the .regions file / GetDP mappings
143 - entities : used for meshing constraints (curve/surface tags)
145 {
146 "dim": <int>, # 0, 1 or 2
147 "phys_tag": <int>, # physical group tag
148 "entities": [int], # list of entity tags (points/curves/surfaces)
149 }
151 The method is idempotent and can be called multiple times after
152 gmsh.model.occ.synchronize().
153 """
155 self.physical_surfaces_by_name.clear()
156 self.physical_curves_by_name.clear()
157 self.physical_points_by_name.clear()
159 phys_groups = gmsh.model.getPhysicalGroups()
161 for dim, pg_tag in phys_groups:
163 name = gmsh.model.getPhysicalName(dim, pg_tag)
164 entity_tags: List[int] = gmsh.model.getEntitiesForPhysicalGroup(dim, pg_tag)
166 meta: Dict[str, Any] = { "dim": dim, "phys_tag": pg_tag, "entities": entity_tags }
167 if dim == 2:
168 self.physical_surfaces_by_name[name] = meta
170 elif dim == 1:
171 self.physical_curves_by_name[name] = meta
173 elif dim == 0:
174 self.physical_points_by_name[name] = meta
178 def _compute_hts_and_layer_mesh_parameters(self) -> None:
179 """
180 Compute the HTS base mesh sizes and approximate vertical
181 discretisation parameters for the stacked metal layers.
183 HTS:
184 - Uses HTS_n_elem_width / HTS_n_elem_thickness to define
185 base element sizes hts_base_size_x / hts_base_size_y.
187 Other layers (Substrate, Silver*, CopperTop/Bottom):
188 - Derive a default element count from the HTS vertical
189 reference size h_ref = hts_base_size_y and the layer
190 thickness (from geometry bounding boxes).
191 - Multiply that default count by the YAML *_elem_scale
192 (substrate_elem_scale, silver_elem_scale,
193 copper_elem_scale) when provided.
194 - Store the final counts and corresponding base sizes in:
195 * self.vertical_elems_by_layer[layer_name]
196 * self.layer_base_size_y[layer_name]
197 """
199 # HTS geometry and mesh parameters
200 s = self.cc_strand
202 W = s.HTS_width
203 T_hts = s.HTS_thickness
204 n_w = self.cc_mesh.HTS_n_elem_width
205 n_t = self.cc_mesh.HTS_n_elem_thickness
207 # Canonical node count along width: used for HTS and for all
208 # tape-internal horizontal boundaries.
209 self.horizontal_nodes_x = int(n_w) + 1
211 # Base HTS element sizes.
212 self.hts_base_size_x = W / float(n_w)
213 self.hts_base_size_y = T_hts / float(n_t)
215 # Initialise per-layer containers
216 self.vertical_elems_by_layer = {}
217 self.layer_base_size_y = {}
219 # Register HTS itself as reference.
220 self.vertical_elems_by_layer["HTS"] = int(n_t)
221 self.layer_base_size_y["HTS"] = self.hts_base_size_y
223 # Reference vertical size from HTS.
224 h_ref = self.hts_base_size_y
226 # Definition of stacked layers and which YAML scale to use.
227 # CopperLeft/Right are not part of the vertical stack and are
228 # deliberately ignored here.
229 layer_definitions = {
230 "Substrate": ("substrate_elem_scale", "substrate layer"),
231 "SilverBottom": ("silver_elem_scale", "silver bottom layer"),
232 "SilverTop": ("silver_elem_scale", "silver top layer"),
233 "CopperBottom": ("copper_elem_scale", "copper bottom layer"),
234 "CopperTop": ("copper_elem_scale", "copper top layer"),
235 }
237 for layer_name, (scale_attr, layer_desc) in layer_definitions.items():
239 meta = self.physical_surfaces_by_name.get(layer_name)
241 # Safety block for optional layers that may be missing.
242 if meta is None:
243 continue
245 surface_tags = list(meta.get("entities", []))
248 # Assumption: layers are axis-aligned (x = width, y = thickness).
249 # We infer layer thickness from the y-extent of OCC bounding boxes.
250 #
251 # Note on *_elem_scale semantics:
252 # *_elem_scale scales the *number of elements* (larger -> finer,
253 # smaller -> coarser), not the characteristic length directly.
254 #
255 # Compute thickness in the stacking (y) direction from the union of
256 # all surfaces in this layer.
258 min_y = float("inf")
259 max_y = float("-inf")
261 for surf_tag in surface_tags:
262 _, ymin, _, _, ymax, _ = gmsh.model.getBoundingBox(2, surf_tag)
263 if ymin < min_y:
264 min_y = ymin
265 if ymax > max_y:
266 max_y = ymax
268 thickness = max_y - min_y
270 # Default element count if we used the same vertical size
271 # as the HTS: n_default ~ thickness / h_ref.
272 n_default = max(1, int(round(thickness / h_ref)))
274 # Apply YAML scale factor on the element count, if provided.
275 scale_value = getattr(self.cc_mesh, scale_attr, None)
276 if scale_value is None or scale_value <= 0.0:
277 n_layer = n_default
278 scale_used = 1.0
280 else:
281 n_layer = max(1, int(round(float(scale_value) * n_default)))
282 scale_used = float(scale_value)
284 # Corresponding base element size in this layer.
285 h_layer = thickness / float(n_layer)
287 # Store results
288 self.vertical_elems_by_layer[layer_name] = n_layer
289 self.layer_base_size_y[layer_name] = h_layer
293 def _apply_air_boundary_mesh_size(self) -> None:
294 """
295 Set the air mesh size from the outer air boundary curve only
296 (physical group 'Air_Outer').
298 The size is:
299 h_air = air_boundary_mesh_size_ratio * HTS_width
301 With MeshSizeExtendFromBoundary=1, this boundary size propagates
302 naturally into the air region.
304 Note:
305 - If a background mesh field is active, it may override point-based
306 sizes in parts of the domain. We still set boundary sizing to keep
307 behaviour robust as a fallback/default when fields are absent.
308 """
310 ratio = getattr(self.cc_mesh, "air_boundary_mesh_size_ratio", None)
311 if ratio is None or float(ratio) <= 0.0:
312 return
314 if getattr(self.cc_strand, "HTS_width", None) is None:
315 raise ValueError("HTS_width is not defined; cannot compute air boundary size.")
317 h_air = float(ratio) * float(self.cc_strand.HTS_width)
319 meta = self.physical_curves_by_name.get("Air_Outer")
320 if meta is None:
321 logger.warning("[Mesh-CC] No 'Air_Outer' physical group found.")
322 return
324 curve_tags = list(meta.get("entities", []))
326 point_tags: list[int] = []
327 for ctag in curve_tags:
328 down_tags, _ = gmsh.model.getAdjacencies(1, ctag) # points on curve
329 point_tags.extend(down_tags)
331 point_tags = list(dict.fromkeys(point_tags))
333 gmsh.model.mesh.setSize([(0, ptag) for ptag in point_tags], h_air)
335 # Make boundary sizing propagate into air
336 gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 1)
337 gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 1)
339 logger.info("[Mesh-CC] Air outer boundary size set: h_air=%.6g (ratio=%.6g)", h_air, float(ratio))
344 def _apply_uniform_copper_mesh_field(self) -> None:
345 """
346 Uniform copper sizing with a single knob, and *no* copper-to-air
347 transition logic.
349 - Copper size: controlled by copper_elem_scale (relative to HTS base size).
350 - Air size: optionally set as a constant target value (h_air) inside the
351 background field, with an additional boundary-based air knob handled
352 separately on 'Air_Outer'.
353 - Background field: Min( constant_air , copper_inside_only )
355 Result:
356 - "Uniform copper" here means the *target size inside copper* is constant.
357 Shared boundaries may still be influenced by transfinite node-count
358 constraints from neighbouring layers.
359 - Air is controlled by a coarse constant target size (or by boundary sizing
360 if the field is not used / as a fallback).
361 """
363 copper_layers = ["CopperTop", "CopperBottom", "CopperLeft", "CopperRight"]
365 copper_surface_tags: list[int] = []
366 for name in copper_layers:
367 copper_surface_tags.extend(self._get_layer_surface_tags(name))
368 copper_surface_tags = list(dict.fromkeys(copper_surface_tags))
370 if not copper_surface_tags:
371 logger.info("[Mesh-CC] No copper surfaces found; skipping copper mesh field.")
372 return
374 # --------------------------------------------------------------
375 # Copper size from single knob: copper_elem_scale
376 # --------------------------------------------------------------
377 scale_value = getattr(self.cc_mesh, "copper_elem_scale", None)
378 if scale_value is None or float(scale_value) <= 0.0:
379 scale_value = 1.0
380 scale_value = float(scale_value)
382 if self.hts_base_size_x is None:
383 raise RuntimeError(
384 "hts_base_size_x not computed; cannot derive copper target size."
385 )
387 h_base = float(self.hts_base_size_x)
388 h_cu = h_base / scale_value
390 # --------------------------------------------------------------
391 # Air size from outer-circle knob: air_boundary_mesh_size_ratio
392 # --------------------------------------------------------------
393 ratio = getattr(self.cc_mesh, "air_boundary_mesh_size_ratio", None)
394 if ratio is not None and getattr(self.cc_strand, "HTS_width", None) is not None:
395 h_air = float(ratio) * float(self.cc_strand.HTS_width)
396 else:
397 # Fallback if knob is not present
398 h_air = 20.0 * h_base
400 # --------------------------------------------------------------
401 # Field A: constant air size everywhere
402 # --------------------------------------------------------------
403 f_air = gmsh.model.mesh.field.add("MathEval")
404 gmsh.model.mesh.field.setString(f_air, "F", f"{h_air}")
406 # --------------------------------------------------------------
407 # Field B: constant size inside copper only (Restrict).
408 # Restrict applies the copper size only inside the listed copper surfaces;
409 # outside those surfaces it returns a very large value, so Min(air, restrict)
410 # naturally selects the air target in the air domain.
411 # --------------------------------------------------------------
412 f_cu_const = gmsh.model.mesh.field.add("MathEval")
413 gmsh.model.mesh.field.setString(f_cu_const, "F", f"{h_cu}")
415 f_cu_inside = gmsh.model.mesh.field.add("Restrict")
416 gmsh.model.mesh.field.setNumber(f_cu_inside, "InField", f_cu_const)
417 gmsh.model.mesh.field.setNumbers(f_cu_inside, "SurfacesList", copper_surface_tags)
419 # --------------------------------------------------------------
420 # Background = min(air everywhere, copper inside)
421 # --------------------------------------------------------------
422 f_min = gmsh.model.mesh.field.add("Min")
423 gmsh.model.mesh.field.setNumbers(f_min, "FieldsList", [f_air, f_cu_inside])
424 gmsh.model.mesh.field.setAsBackgroundMesh(f_min)
426 # Let fine boundary sizes propagate into neighbouring domains (air grading).
427 gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 1)
428 gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 1)
430 logger.info(
431 "[Mesh-CC] Air+copper constant background mesh: "
432 "h_air=%.6g (ratio=%s), h_cu=%.6g (scale=%.6g, h_base=%.6g)",
433 h_air,
434 str(ratio),
435 h_cu,
436 scale_value,
437 h_base,
438 )
447 # ------------------------------------------------------------------ #
448 # Internal helpers: access to physical surfaces
449 # ------------------------------------------------------------------ #
451 def _get_layer_surface_tags(self, layer_name: str) -> List[int]:
452 """
453 Return surface tags for a given 2D layer physical group.
455 Parameters
456 ----------
457 layer_name : str
458 Name of the 2D physical group, e.g. "Substrate", "SilverTop".
460 Returns
461 -------
462 list of int
463 List of surface tags
464 """
465 meta = self.physical_surfaces_by_name.get(layer_name)
466 if not meta:
467 return []
468 return list(meta.get("entities", []))
472 def _get_hts_surface_tags(self) -> List[int]:
473 """
474 Return surface tags for HTS-related regions.
476 We treat:
477 - "HTS" : all HTS filaments together (monolithic or striated).
478 - "HTS_Grooves" : groove surfaces (if present, created as
479 "HTS_Grooves" physical group by geometry).
481 as a single logical HTS layer for transfinite meshing. This
482 ensures that filaments and grooves (if any) share a consistent
483 structured discretisation.
484 """
485 tags: List[int] = []
487 meta_hts = self.physical_surfaces_by_name.get("HTS")
488 if meta_hts:
489 tags.extend(meta_hts.get("entities", []))
491 meta_grooves = self.physical_surfaces_by_name.get("HTS_Grooves")
492 if meta_grooves:
493 tags.extend(meta_grooves.get("entities", []))
495 # De-duplicate while preserving order.
496 return list(dict.fromkeys(tags))
500 def _apply_transfinite_to_hts(self) -> None:
501 """
502 Build a structured (transfinite) mesh for HTS filaments, and keep
503 groove regions (HTS_Grooves) unstructured.
505 - Filaments ("HTS" physical group):
506 * Horizontal edges: transfinite with 'Bump' distribution.
507 * Vertical edges: uniform spacing.
508 * Surfaces: transfinite + recombine -> quadrilateral mesh.
510 - Grooves ("HTS_Grooves" physical group):
511 * Boundary curves still get transfinite node counts so that
512 they align with neighbouring HTS and metal layers.
513 * BUT we do NOT call setTransfiniteSurface / setRecombine,
514 so the interior of grooves is meshed as unstructured triangles.
515 """
516 # All HTS-related 2D surfaces (filaments + grooves).
517 surface_tags = self._get_hts_surface_tags()
519 # ------------------------------------------------------------------
520 # Step 1: total filament width and reference Δx.
521 # ------------------------------------------------------------------
522 filament_meta = self.physical_surfaces_by_name.get("HTS", {})
523 filament_surface_tags: List[int] = list(
524 filament_meta.get("entities", [])
525 )
527 equiv_width_filaments = 0.0
528 for surf_tag in filament_surface_tags:
529 xmin, _, _, xmax, _, _ = gmsh.model.getBoundingBox(2, surf_tag)
530 dx = xmax - xmin
531 if dx > 0.0:
532 equiv_width_filaments += dx
534 n_tot = float(self.cc_mesh.HTS_n_elem_width)
536 # Reset per-surface maps by default.
537 self.hts_nodes_x_filaments = {}
538 self.hts_nodes_x_grooves = {}
540 if equiv_width_filaments > 0.0 and n_tot > 0.0:
541 # Reference horizontal cell size in superconducting material:
542 # Δx_ref = W_SC / N_tot
543 # where W_SC is the sum of all filament widths.
544 h_ref_x = equiv_width_filaments / n_tot
545 self.equiv_hts_filament_width = equiv_width_filaments
546 self.hts_ref_size_x_filaments = h_ref_x
548 logger.debug(
549 "[Mesh-CC] HTS equivalent filament width = %.6g, "
550 "HTS reference horizontal size (filaments) = %.6g",
551 equiv_width_filaments,
552 h_ref_x,
553 )
555 # ------------------------------------------------------------------
556 # Step 2: horizontal element densities in filaments and grooves.
557 # ------------------------------------------------------------------
558 # Elements per unit length in filaments:
559 # λ_fil = N_tot / W_SC (≈ 1 / Δx_ref)
560 lambda_fil = n_tot / equiv_width_filaments
562 # Groove density ratio (fewer elements in grooves).
563 groove_ratio = getattr(
564 self.cc_mesh, "HTS_groove_elem_density_ratio", None
565 )
566 if groove_ratio is None or groove_ratio <= 0.0:
567 groove_ratio = 0.25
569 lambda_groove = groove_ratio * lambda_fil
571 self.hts_density_x_filaments = lambda_fil
572 self.hts_density_x_grooves = lambda_groove
574 logger.debug(
575 "[Mesh-CC] HTS filament density (lambda_fil) = %.6g, "
576 "groove density ratio = %.6g, "
577 "groove density (lambda_groove) = %.6g",
578 lambda_fil,
579 groove_ratio,
580 lambda_groove,
581 )
583 # ------------------------------------------------------------------
584 # Step 3: per-surface horizontal node counts.
585 # ------------------------------------------------------------------
586 groove_min_elems = getattr(
587 self.cc_mesh, "HTS_groove_min_elems", None
588 )
589 if groove_min_elems is None or groove_min_elems < 1:
590 groove_min_elems = 1
592 # --- Filaments: N_i ≈ λ_fil * width_i ---
593 for surf_tag in filament_surface_tags:
594 xmin, _, _, xmax, _, _ = gmsh.model.getBoundingBox(2, surf_tag)
595 width = xmax - xmin
596 if width <= 0.0:
597 continue
599 n_elem_f = max(1, int(round(lambda_fil * width)))
600 n_nodes_x_f = n_elem_f + 1 # nodes = elements + 1
602 self.hts_nodes_x_filaments[surf_tag] = n_nodes_x_f
604 logger.debug(
605 "[Mesh-CC] HTS filament surface %d: width = %.6g, "
606 "n_elem = %d, n_nodes_x = %d",
607 surf_tag,
608 width,
609 n_elem_f,
610 n_nodes_x_f,
611 )
613 # --- Grooves: fewer elements per unit length ---
614 groove_meta = self.physical_surfaces_by_name.get(
615 "HTS_Grooves", {}
616 )
617 groove_surface_tags: List[int] = list(
618 groove_meta.get("entities", [])
619 )
621 for surf_tag in groove_surface_tags:
622 xmin, _, _, xmax, _, _ = gmsh.model.getBoundingBox(2, surf_tag)
623 width = xmax - xmin
624 if width <= 0.0:
625 continue
627 n_elem_g = int(round(lambda_groove * width))
628 n_elem_g = max(groove_min_elems, n_elem_g)
629 n_nodes_x_g = n_elem_g + 1
631 self.hts_nodes_x_grooves[surf_tag] = n_nodes_x_g
633 logger.debug(
634 "[Mesh-CC] HTS groove surface %d: width = %.6g, "
635 "n_elem = %d, n_nodes_x = %d",
636 surf_tag,
637 width,
638 n_elem_g,
639 n_nodes_x_g,
640 )
642 else:
643 # Fallback: no HTS-specific filament logic.
644 self.equiv_hts_filament_width = None
645 self.hts_ref_size_x_filaments = None
646 self.hts_density_x_filaments = None
647 self.hts_density_x_grooves = None
648 self.hts_nodes_x_filaments = {}
649 self.hts_nodes_x_grooves = {}
651 # ------------------------------------------------------------------
652 # Step 4: Apply transfinite curves.
653 #
654 # - Horizontal curves (dx >= dy): 'Bump' with surf_n_nodes_x nodes.
655 # - Vertical curves (dy > dx): uniform 'Progression' with n_nodes_y.
656 #
657 # Then:
658 # - setTransfiniteSurface + Recombine ONLY for filament surfaces.
659 # - Grooves: NO setTransfiniteSurface / NO Recombine -> unstructured.
660 # ------------------------------------------------------------------
662 # Global vertical node count (same for all HTS surfaces).
663 n_t = self.cc_mesh.HTS_n_elem_thickness
664 n_nodes_y = int(n_t) + 1
666 # Fallback horizontal node count if per-surface info is missing.
667 n_w = self.cc_mesh.HTS_n_elem_width
668 fallback_n_nodes_x = (
669 self.horizontal_nodes_x
670 if self.horizontal_nodes_x is not None
671 else int(n_w) + 1
672 )
674 # Bump coefficient for clustering along width.
675 bump_coef = getattr(self.cc_mesh, "bump_coef", None)
676 if bump_coef is None or bump_coef <= 0.0:
677 bump_coef = 0.01
679 # Re-build sets of filament and groove surfaces for classification.
680 filament_meta = self.physical_surfaces_by_name.get("HTS", {})
681 filament_surface_set = set(filament_meta.get("entities", []))
683 groove_meta = self.physical_surfaces_by_name.get("HTS_Grooves", {})
684 groove_surface_set = set(groove_meta.get("entities", []))
686 # Track all HTS boundary curves so other routines can skip them.
687 self.hts_boundary_curves = set()
689 for surf_tag in surface_tags:
690 # Horizontal node count for this specific surface.
691 if surf_tag in self.hts_nodes_x_filaments:
692 surf_n_nodes_x = self.hts_nodes_x_filaments[surf_tag]
693 elif surf_tag in self.hts_nodes_x_grooves:
694 surf_n_nodes_x = self.hts_nodes_x_grooves[surf_tag]
695 else:
696 surf_n_nodes_x = fallback_n_nodes_x
698 # Get boundary curves of this HTS region.
699 boundary = gmsh.model.getBoundary(
700 [(2, surf_tag)], oriented=False, recursive=False
701 )
702 curve_tags = [c[1] for c in boundary if c[0] == 1]
703 curve_tags = list(dict.fromkeys(curve_tags))
705 for ctag in curve_tags:
706 xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(1, ctag)
707 dx = abs(xmax - xmin)
708 dy = abs(ymax - ymin)
710 # Mark as HTS boundary curve.
711 self.hts_boundary_curves.add(ctag)
713 if dx >= dy:
714 # Horizontal curve -> controls columns (along tape width).
715 gmsh.model.mesh.setTransfiniteCurve(
716 ctag,
717 surf_n_nodes_x,
718 "Bump",
719 float(bump_coef),
720 )
721 else:
722 # Vertical curve -> controls rows (through thickness).
723 gmsh.model.mesh.setTransfiniteCurve(
724 ctag,
725 n_nodes_y,
726 "Progression",
727 1.0,
728 )
730 # --- Surface-level transfinite / recombine behaviour ---
731 if surf_tag in filament_surface_set:
732 # Filaments: fully structured quads.
733 gmsh.model.mesh.setTransfiniteSurface(surf_tag)
734 gmsh.model.mesh.setRecombine(2, surf_tag)
735 elif surf_tag in groove_surface_set:
736 # Grooves: DO NOT set transfinite surface and DO NOT recombine.
737 # The 2D mesher will generate an unstructured triangular mesh
738 # inside, but boundary node distributions are still honoured.
739 continue
740 else:
741 # Should not happen (all HTS surfaces are either filament
742 # or groove), but keep a safe fallback.
743 gmsh.model.mesh.setTransfiniteSurface(surf_tag)
746 def _apply_horizontal_boundary_subdivisions_for_side_copper(self) -> None:
747 """
748 Apply transfinite subdivisions along the tape width direction for
749 CopperLeft / CopperRight so that their horizontal boundaries use the
750 same column pattern as the HTS / stacked layers.
751 """
753 # Canonical horizontal node count along the tape width.
754 n_nodes_horizontal = self.horizontal_nodes_x
755 if n_nodes_horizontal is None or n_nodes_horizontal <= 1:
756 # Fallback to HTS_n_elem_width + 1 if something went wrong.
757 n_w = self.cc_mesh.HTS_n_elem_width
758 n_nodes_horizontal = int(n_w) + 1
760 # Air surfaces: used to detect copper-air edges that we want
761 # to keep coarser than the internal tape interfaces.
762 air_surface_tags = set(self._get_layer_surface_tags("Air"))
764 layer_names = ["CopperLeft", "CopperRight"]
766 for layer_name in layer_names:
767 surface_tags = self._get_layer_surface_tags(layer_name)
769 for surf_tag in surface_tags:
770 boundary = gmsh.model.getBoundary(
771 [(2, surf_tag)], oriented=False, recursive=False
772 )
773 curve_tags = [c[1] for c in boundary if c[0] == 1]
774 curve_tags = list(dict.fromkeys(curve_tags))
776 for ctag in curve_tags:
777 xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(1, ctag)
778 dx = abs(xmax - xmin)
779 dy = abs(ymax - ymin)
781 # Only curves that run mainly along x (tape width).
782 if dx >= dy:
783 # Check if this horizontal edge touches air.
784 higher_dim_tags, _ = gmsh.model.getAdjacencies(1, ctag)
786 if any(s in air_surface_tags for s in higher_dim_tags):
787 # Copper–air edge: keep this much coarser than internal
788 # tape interfaces to avoid over-refining the surrounding air.
789 # Heuristic: ~2.5% of the internal column count.
790 n_outer = max(
791 2, int(round(0.025 * n_nodes_horizontal))
792 )
793 gmsh.model.mesh.setTransfiniteCurve(
794 ctag,
795 n_outer,
796 "Progression",
797 1.0,
798 )
799 continue
801 # Internal horizontal interface: align with tape columns.
802 gmsh.model.mesh.setTransfiniteCurve(
803 ctag,
804 n_nodes_horizontal,
805 "Progression",
806 1.0,
807 )
810 def _apply_vertical_boundary_subdivisions_for_side_copper(self) -> None:
811 """
812 Apply transfinite node counts to vertical boundary curves of
813 CopperLeft / CopperRight that are adjacent to the Air region.
814 """
815 # Requires _compute_hts_and_layer_mesh_parameters() to have run
816 # (hts_base_size_y set).
817 h_ref = self.hts_base_size_y
819 # YAML scale for copper; <= 0 or None -> use 1.0 (no extra scaling)
820 scale_value = getattr(self.cc_mesh, "copper_elem_scale", None)
821 if scale_value is None or scale_value <= 0.0:
822 scale_value = 1.0
824 # Air surfaces, to detect copper-air vertical edges
825 air_surface_tags = set(self._get_layer_surface_tags("Air"))
827 layer_names = ["CopperLeft", "CopperRight"]
829 for layer_name in layer_names:
830 surface_tags = self._get_layer_surface_tags(layer_name)
832 for surf_tag in surface_tags:
833 boundary = gmsh.model.getBoundary( [(2, surf_tag)], oriented=False, recursive=False )
834 curve_tags = [ c[1] for c in boundary if c[0] == 1 ]
835 curve_tags = list(dict.fromkeys(curve_tags))
837 for ctag in curve_tags:
838 xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(1, ctag)
839 dx = abs(xmax - xmin)
840 dy = abs(ymax - ymin)
842 # We only care about vertical curves here.
843 if dy <= dx:
844 continue
846 # Check if this vertical curve touches Air.
847 higher_dim_tags, _ = gmsh.model.getAdjacencies(1, ctag)
848 if not any(s in air_surface_tags for s in higher_dim_tags):
849 # Internal copper-tape edge: stacked-layer logic handles it.
850 continue
852 thickness = dy
853 n_default = max(1, int(round(thickness / h_ref)))
854 n_layer = max(1, int(round(scale_value * n_default)))
855 n_nodes_vertical = n_layer + 1
857 gmsh.model.mesh.setTransfiniteCurve(
858 ctag,
859 n_nodes_vertical,
860 "Progression",
861 1.0,
862 )
865 def _apply_vertical_boundary_subdivisions_to_non_hts_layers(self) -> None:
866 """
867 Apply transfinite node counts to boundary curves of non-HTS, non-copper
868 stacked layers (Substrate, Silver*), based on:
869 - vertical_elems_by_layer[layer] for vertical curves (thickness dir),
870 - self.horizontal_nodes_x for horizontal curves (width dir).
872 Copper is intentionally excluded here to preserve more uniform copper
873 node spacing controlled by the background sizing field.
874 """
876 # Horizontal node count along the tape width.
877 n_nodes_horizontal = self.horizontal_nodes_x
879 # Set of all HTS-related surfaces: used to avoid overriding
880 # transfinite settings on HTS boundary curves.
881 hts_surface_tags = set(self._get_hts_surface_tags())
883 # Non-HTS stacked layers that we want to control in the vertical
884 # (thickness) direction and in the horizontal (width) direction
885 # via the canonical column count.
886 layer_names = [
887 "Substrate",
888 "SilverBottom",
889 "SilverTop",
890 ]
893 for layer_name in layer_names:
894 n_layer = self.vertical_elems_by_layer.get(layer_name)
896 if n_layer is None:
897 continue
899 surface_tags = self._get_layer_surface_tags(layer_name)
901 # Number of nodes on vertical curves = number of elements + 1
902 n_nodes_vertical = int(n_layer) + 1
904 for surf_tag in surface_tags:
905 # Get boundary curves of this surface
906 boundary = gmsh.model.getBoundary( [(2, surf_tag)], oriented=False, recursive=False )
907 curve_tags = [c[1] for c in boundary if c[0] == 1]
908 curve_tags = list(dict.fromkeys(curve_tags))
910 for ctag in curve_tags:
911 # If this curve is shared with an HTS surface,
912 # skip it so the HTS progression is preserved.
913 higher_dim_tags, _ = gmsh.model.getAdjacencies(1, ctag)
914 if any(s in hts_surface_tags for s in higher_dim_tags):
915 logger.debug(
916 "[Mesh-CC]'%s' curve %d touches HTS; "
917 "skipping non-HTS transfinite constraint.",
918 layer_name,
919 ctag,
920 )
921 continue
923 xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(1, ctag)
924 dx = abs(xmax - xmin)
925 dy = abs(ymax - ymin)
927 if dy > dx:
928 # Predominantly vertical curve: control thickness direction
929 gmsh.model.mesh.setTransfiniteCurve(
930 ctag, n_nodes_vertical, "Progression", 1.0
931 )
932 else:
933 # Predominantly horizontal curve: enforce common columns.
934 # For SilverTop outer boundary, use a Bump distribution to mimic HTS clustering.
935 if layer_name == "SilverTop":
936 bump_coef = getattr(self.cc_mesh, "bump_coef", None)
938 if bump_coef is None or float(bump_coef) <= 0.0:
939 # Fallback: reuse HTS bump coefficient (or default)
940 bump_coef = getattr(self.cc_mesh, "bump_coef", None)
942 if bump_coef is None or float(bump_coef) <= 0.0:
943 bump_coef = 0.01
945 gmsh.model.mesh.setTransfiniteCurve(
946 ctag, n_nodes_horizontal, "Bump", float(bump_coef)
947 )
948 else:
949 gmsh.model.mesh.setTransfiniteCurve(
950 ctag, n_nodes_horizontal, "Progression", 1.0
951 )
956 def _apply_progression_on_substrate_vertical_sides(self) -> None:
957 """
958 Apply a 'Progression' distribution on the vertical side boundary
959 curves of the substrate, clustering nodes toward the HTS side (top in y).
960 """
961 r = getattr(self.cc_mesh, "substrate_side_progression", None)
962 if r is None:
963 return
965 r = float(r)
966 if r <= 1.0:
967 return
969 n_layer = self.vertical_elems_by_layer.get("Substrate")
970 if n_layer is None:
971 return
972 n_nodes_vertical = int(n_layer) + 1
974 substrate_surface_tags = self._get_layer_surface_tags("Substrate")
975 if not substrate_surface_tags:
976 return
978 # All HTS-related surfaces (HTS + grooves) to detect the top interface curve.
979 hts_surface_tags = set(self._get_hts_surface_tags())
981 # Collect candidate curves: substrate boundary curves that are vertical
982 # and do NOT touch HTS (so we keep the HTS/substrate interface unchanged).
983 candidate_curve_tags: list[int] = []
985 for surf_tag in substrate_surface_tags:
986 boundary = gmsh.model.getBoundary(
987 [(2, surf_tag)], oriented=False, recursive=False
988 )
989 curve_tags = [c[1] for c in boundary if c[0] == 1]
990 curve_tags = list(dict.fromkeys(curve_tags))
992 for ctag in curve_tags:
993 # Skip curves that touch HTS surfaces (top interface).
994 higher_dim_tags, _ = gmsh.model.getAdjacencies(1, ctag)
995 if any(s in hts_surface_tags for s in higher_dim_tags):
996 continue
998 xmin, ymin, _, xmax, ymax, _ = gmsh.model.getBoundingBox(1, ctag)
999 dx = abs(xmax - xmin)
1000 dy = abs(ymax - ymin)
1002 # Only vertical-ish curves (substrate left/right sides).
1003 if dy <= dx:
1004 continue
1006 candidate_curve_tags.append(ctag)
1008 candidate_curve_tags = list(dict.fromkeys(candidate_curve_tags))
1009 if not candidate_curve_tags:
1010 logger.info(
1011 "[Mesh-CC] No substrate vertical side curves found for Progression."
1012 )
1013 return
1015 # Apply progression with refinement toward *top end* (larger y).
1016 for ctag in candidate_curve_tags:
1017 endpoints = gmsh.model.getBoundary(
1018 [(1, ctag)], oriented=True, recursive=False
1019 )
1020 point_tags = [p[1] for p in endpoints if p[0] == 0]
1022 if len(point_tags) < 2:
1023 continue
1025 p_start = point_tags[0]
1026 p_end = point_tags[-1]
1028 _, y0, _ = gmsh.model.getValue(0, p_start, [])
1029 _, y1, _ = gmsh.model.getValue(0, p_end, [])
1031 # Progression behaviour:
1032 # - ratio > 1 -> smallest segments near start
1033 # - ratio < 1 -> smallest segments near end
1034 #
1035 # Note: refinement direction depends on curve orientation.
1036 # We infer which endpoint is "top" by comparing point y-coordinates
1037 # and invert the progression ratio accordingly.
1038 #
1039 # We want smallest segments near the *top* (higher y).
1040 if y0 >= y1:
1041 prog = r
1042 else:
1043 prog = 1.0 / r
1045 gmsh.model.mesh.setTransfiniteCurve(
1046 ctag,
1047 n_nodes_vertical,
1048 "Progression",
1049 float(prog),
1050 )
1055 # ------------------------------------------------------------------ #
1056 # Public API
1057 # ------------------------------------------------------------------ #
1059 def generate_mesh(self, geom_folder: str) -> None:
1061 geo_unrolled_path = os.path.join(geom_folder, f"{self.magnet_name}.xao")
1063 # The .xao file is produced by GeometryConductorAC_CC and must contain the
1064 # expected physical group names used throughout this meshing routine
1065 # (HTS, optional HTS_Grooves, metal layers, Air, Air_Outer, etc.).
1066 gmsh.open(geo_unrolled_path)
1067 occ.synchronize()
1069 logger.debug("[Mesh-CC] OCC model synchronized after loading geometry file.")
1071 # Collect all physical groups into helper dictionaries that will
1072 # be reused in later steps (mesh fields, regions, ...).
1073 self._collect_physical_groups()
1075 # HTS base mesh sizes and vertical element counts.
1076 self._compute_hts_and_layer_mesh_parameters()
1078 # Transfinite meshing for HTS (filaments + grooves or monolithic).
1079 self._apply_transfinite_to_hts()
1081 # Apply transfinite constraints only to non-copper stacked layers
1082 # (Substrate, Silver*) so we do not distort copper node spacing.
1083 self._apply_vertical_boundary_subdivisions_to_non_hts_layers()
1085 # Refine substrate vertical sides near the HTS side using Progression.
1086 # This overrides the uniform 'Progression, 1.0' set above for those curves.
1087 self._apply_progression_on_substrate_vertical_sides()
1090 # Copper is handled by a uniform size field (no transfinite on copper curves).
1091 self._apply_uniform_copper_mesh_field()
1095 # -----------------------------------------------------------------
1096 # Global mesh size scaling and actual 2D mesh generation
1097 # -----------------------------------------------------------------
1098 scaling_global = self.cc_mesh.scaling_global
1100 # MeshSizeFactor scales characteristic lengths (point sizes, background
1101 # field values, etc.). Transfinite constraints are node *counts* and are
1102 # not affected by this factor.
1103 gmsh.option.setNumber("Mesh.MeshSizeFactor", float(scaling_global))
1105 # Boundary-based air sizing knob (fallback/default behaviour).
1106 # Note: a background field (if active) can override sizes in parts of the
1107 # domain; this still provides a robust air sizing rule when fields are not used.
1108 self._apply_air_boundary_mesh_size()
1110 # -----------------------------------------------------------------
1111 # Generate 1D first to "freeze" transfinite curve discretisations,
1112 # then generate 2D. This avoids Gmsh re-meshing/compacting 1D curves
1113 # when strong background fields are present.
1114 # -----------------------------------------------------------------
1115 gmsh.model.mesh.generate(1)
1116 gmsh.model.mesh.generate(2)
1118 # Generate cohomology basis functions (cuts) for transport current handling with h-phi-formulation
1119 self.generate_cuts()
1121 # Generate region file for easy access to physical groups from, e.g., the .pro template
1122 self.generate_regions_file()
1124 # Save the mesh to <magnet_name>.msh and optionally launch the GUI.
1125 gui = getattr(self.fdm.run, "launch_gui", False)
1126 self.save_mesh(gui=gui)
1128 logger.info( "[Mesh-CC] Mesh generation and saving complete. Mesh file: '%s'.", self.mesh_file )
1131 # -----------------------------------------------------------------
1132 # Generation of .regions file
1133 # -----------------------------------------------------------------
1135 def generate_regions_file(self) -> None:
1136 """
1137 Construct the .regions YAML file for the 2D CACCC conductor model.
1138 """
1139 # Initialise the RegionsModel container
1140 rm = RegionsModel()
1142 # Powered - HTS
1143 rm.powered["HTS"] = RegionsModelFiQuS.Powered()
1144 powered = rm.powered["HTS"]
1146 # Initialise list-valued fields
1147 for fld in ("vol", "surf", "surf_th", "surf_in", "surf_out", "curve", "surf_insul"):
1148 getattr(powered, fld).names = []
1149 getattr(powered, fld).numbers = []
1151 # Terminal volume placeholders (not used in 2D but required structurally)
1152 powered.vol_in = RegionsModelFiQuS.Region(name=None, number=None)
1153 powered.vol_out = RegionsModelFiQuS.Region(name=None, number=None)
1155 # Convention: for 2D models we store surface physical tags in the "vol"
1156 # fields to keep the RegionsModel structure compatible with 3D templates
1157 # and existing .pro logic.
1158 if "HTS" in self.physical_surfaces_by_name:
1159 meta = self.physical_surfaces_by_name["HTS"]
1160 powered.vol.names = ["HTS"]
1161 powered.vol.numbers = [meta["phys_tag"]]
1163 # HTS boundary curves
1164 for edge_name in ["HTS_Upper", "HTS_Lower", "HTS_LeftEdge", "HTS_RightEdge"]:
1165 if edge_name in self.physical_curves_by_name:
1166 meta = self.physical_curves_by_name[edge_name]
1167 powered.surf_th.names.append(edge_name)
1168 powered.surf_th.numbers.append(meta["phys_tag"])
1171 # Induced - Substrate + Stabilizer layers
1172 rm.induced["Stabilizer"] = RegionsModelFiQuS.Induced()
1173 induced = rm.induced["Stabilizer"]
1175 # Initialise all induced list-valued fields
1176 for fld in ("vol", "surf_th", "surf_in", "surf_out", "cochain"):
1177 getattr(induced, fld).names = []
1178 getattr(induced, fld).numbers = []
1180 # Metallic layer volumes (2D surfaces)
1181 metal_layers = [
1182 "Substrate",
1183 "SilverTop", "SilverBottom",
1184 "CopperTop", "CopperBottom",
1185 "CopperLeft", "CopperRight",
1186 ]
1188 for layer in metal_layers:
1189 if layer in self.physical_surfaces_by_name:
1190 meta = self.physical_surfaces_by_name[layer]
1191 induced.vol.names.append(layer)
1192 induced.vol.numbers.append(meta["phys_tag"])
1194 # Outer air boundary shared with the stabilizer
1195 if "Air_Outer" in self.physical_curves_by_name:
1196 meta = self.physical_curves_by_name["Air_Outer"]
1197 induced.surf_out.names = ["Air_Outer"]
1198 induced.surf_out.numbers = [meta["phys_tag"]]
1200 # Metallic layer boundary edges
1201 edge_suffixes = ["Upper", "Lower", "LeftEdge", "RightEdge"]
1203 for layer in metal_layers:
1204 for suffix in edge_suffixes:
1205 name = f"{layer}_{suffix}"
1206 if name in self.physical_curves_by_name:
1207 meta = self.physical_curves_by_name[name]
1208 induced.surf_th.names.append(name)
1209 induced.surf_th.numbers.append(meta["phys_tag"])
1211 # Explicit interface curves between layers
1212 interface_groups = [
1213 "Buffer Layer",
1214 "HTS_Silver_Interface",
1215 "HTS_CopperTop_Interface",
1216 "HTS_CopperLeft_Interface",
1217 "HTS_CopperRight_Interface",
1218 "SilverTop_CopperTop_Interface",
1219 "SilverBottom_CopperBottom_Interface",
1220 "Substrate_SilverBottom_Interface",
1221 "Substrate_CopperBottom_Interface",
1222 "Substrate_CopperLeft_Interface",
1223 "Substrate_CopperRight_Interface",
1224 ]
1226 for name in interface_groups:
1227 if name in self.physical_curves_by_name:
1228 meta = self.physical_curves_by_name[name]
1229 induced.surf_th.names.append(name)
1230 induced.surf_th.numbers.append(meta["phys_tag"])
1233 # Air region
1234 if "Air" in self.physical_surfaces_by_name:
1235 meta = self.physical_surfaces_by_name["Air"]
1236 rm.air.vol.name = "Air"
1237 rm.air.vol.number = meta["phys_tag"]
1239 if "Air_Outer" in self.physical_curves_by_name:
1240 meta = self.physical_curves_by_name["Air_Outer"]
1241 rm.air.surf.name = "Air_Outer"
1242 rm.air.surf.number = meta["phys_tag"]
1244 if "Air_Inner" in self.physical_curves_by_name:
1245 meta = self.physical_curves_by_name["Air_Inner"]
1246 rm.air.line.name = "Air_Inner"
1247 rm.air.line.number = meta["phys_tag"]
1249 if "Gauging_point" in self.physical_points_by_name:
1250 meta = self.physical_points_by_name["Gauging_point"]
1251 rm.air.point.names = ["Gauging_point"]
1252 rm.air.point.numbers = [meta["phys_tag"]]
1254 # Cut for transport current
1255 if "Cut" in self.physical_curves_by_name:
1256 meta = self.physical_curves_by_name["Cut"]
1257 rm.air.cochain.names = ["Cut"]
1258 rm.air.cochain.numbers = [meta["phys_tag"]]
1261 # Write the regions file and log
1262 FilesAndFolders.write_data_to_yaml(self.regions_file, rm.model_dump())
1263 logger.info(f"[Mesh-CC] Regions file written to: {self.regions_file}")
1266 def generate_cuts(self) -> None:
1267 """
1268 Generate cohomology basis functions for the CACCC model.
1269 """
1270 # Transport current cut
1271 air_inner = self.physical_curves_by_name["Air_Inner"]
1272 air_surf = self.physical_surfaces_by_name["Air"]
1274 # Request one 1D cut cycle (homology) and the corresponding 1D cohomology
1275 # basis inside the Air surface. This supports transport-current handling in
1276 # h-phi formulations.
1277 gmsh.model.mesh.addHomologyRequest(
1278 "Homology", domainTags=[air_inner["phys_tag"]], dims=[1]
1279 )
1280 gmsh.model.mesh.addHomologyRequest(
1281 "Cohomology", domainTags=[air_surf["phys_tag"]], dims=[1]
1282 )
1284 # computeHomology() returns tags of generated chains/physical groups.
1285 # The indexing below follows Gmsh's returned ordering for our request pair.
1286 cuts = gmsh.model.mesh.computeHomology()
1288 # Rearrange so that positive currents have deterministic orientation.
1289 gmsh.model.mesh.clearHomologyRequests()
1290 gmsh.plugin.setString(
1291 "HomologyPostProcessing",
1292 "PhysicalGroupsOfOperatedChains",
1293 str(cuts[1][1]),
1294 )
1295 gmsh.plugin.setString(
1296 "HomologyPostProcessing",
1297 "PhysicalGroupsOfOperatedChains2",
1298 str(cuts[0][1]),
1299 )
1300 gmsh.plugin.run("HomologyPostProcessing")
1302 # Store the physical group tag for later use (e.g. writing to the regions file).
1303 # The "+1" matches the physical group id created by the post-processing plugin.
1304 meta: Dict[str, Any] = {"dim": 1, "phys_tag": cuts[1][1] + 1}
1305 self.physical_curves_by_name["Cut"] = meta