add tunable corner-angle splitting for bspline meshing

Author: mdmaas

nbs/palace_dc.ipynb

@@ -345,6 +345,7 @@ def _build_mesh_config(
         curve_fit_layers: list[str] | None = None,
         curve_fit_tolerance_um: float | None = None,
         curve_fit_min_points: int | None = None,
+        curve_fit_corner_angle_deg: float | None = None,
         high_order_elements: bool | None = None,
         high_order_order: int | None = None,
         high_order_optimize: bool | None = None,
@@ -427,6 +428,9 @@ def _build_mesh_config(
             mesh_config.curve_fit_layers = list(existing_config.curve_fit_layers)
             mesh_config.curve_fit_tolerance_um = existing_config.curve_fit_tolerance_um
             mesh_config.curve_fit_min_points = existing_config.curve_fit_min_points
+            mesh_config.curve_fit_corner_angle_deg = (
+                existing_config.curve_fit_corner_angle_deg
+            )
             mesh_config.high_order_elements = existing_config.high_order_elements
             mesh_config.high_order_order = existing_config.high_order_order
             mesh_config.high_order_optimize = existing_config.high_order_optimize
@@ -462,6 +466,8 @@ def _build_mesh_config(
             mesh_config.curve_fit_tolerance_um = curve_fit_tolerance_um
         if curve_fit_min_points is not None:
             mesh_config.curve_fit_min_points = curve_fit_min_points
+        if curve_fit_corner_angle_deg is not None:
+            mesh_config.curve_fit_corner_angle_deg = curve_fit_corner_angle_deg
         if high_order_elements is not None:
             mesh_config.high_order_elements = high_order_elements
         if high_order_order is not None:
@@ -930,6 +936,7 @@ def _generate_mesh_internal(
             curve_fit_layers=mesh_config.curve_fit_layers,
             curve_fit_tolerance_um=mesh_config.curve_fit_tolerance_um,
             curve_fit_min_points=mesh_config.curve_fit_min_points,
+            curve_fit_corner_angle_deg=mesh_config.curve_fit_corner_angle_deg,
             high_order_elements=mesh_config.high_order_elements,
             high_order_order=mesh_config.high_order_order,
             high_order_optimize=mesh_config.high_order_optimize,
@@ -979,6 +986,7 @@ def preview(
         curve_fit_layers: list[str] | None = None,
         curve_fit_tolerance_um: float | None = None,
         curve_fit_min_points: int | None = None,
+        curve_fit_corner_angle_deg: float | None = None,
         high_order_elements: bool | None = None,
         high_order_order: int | None = None,
         high_order_optimize: bool | None = None,
@@ -1006,6 +1014,8 @@ def preview(
             curve_fit_layers: Layer names where spline fitting is applied.
             curve_fit_tolerance_um: Point merge tolerance before fitting.
             curve_fit_min_points: Min contour points required to fit curves.
+            curve_fit_corner_angle_deg: Turn-angle threshold used to split
+                sharp corners from smooth curve-fit segments.
             high_order_elements: Enable high-order geometric mesh elements.
             high_order_order: Polynomial order for high-order elements.
             high_order_optimize: Run gmsh high-order optimization after meshing.
@@ -1040,6 +1050,7 @@ def preview(
             curve_fit_layers=curve_fit_layers,
             curve_fit_tolerance_um=curve_fit_tolerance_um,
             curve_fit_min_points=curve_fit_min_points,
+            curve_fit_corner_angle_deg=curve_fit_corner_angle_deg,
             high_order_elements=high_order_elements,
             high_order_order=high_order_order,
             high_order_optimize=high_order_optimize,
@@ -1076,6 +1087,7 @@ def preview(
                 curve_fit_layers=mesh_config.curve_fit_layers,
                 curve_fit_tolerance_um=mesh_config.curve_fit_tolerance_um,
                 curve_fit_min_points=mesh_config.curve_fit_min_points,
+                curve_fit_corner_angle_deg=mesh_config.curve_fit_corner_angle_deg,
                 high_order_elements=mesh_config.high_order_elements,
                 high_order_order=mesh_config.high_order_order,
                 high_order_optimize=mesh_config.high_order_optimize,
@@ -1106,6 +1118,7 @@ def mesh(
         curve_fit_layers: list[str] | None = None,
         curve_fit_tolerance_um: float | None = None,
         curve_fit_min_points: int | None = None,
+        curve_fit_corner_angle_deg: float | None = None,
         high_order_elements: bool | None = None,
         high_order_order: int | None = None,
         high_order_optimize: bool | None = None,
@@ -1139,6 +1152,8 @@ def mesh(
             curve_fit_layers: Layer names where spline fitting is applied.
             curve_fit_tolerance_um: Point merge tolerance before fitting.
             curve_fit_min_points: Min contour points required to fit curves.
+            curve_fit_corner_angle_deg: Turn-angle threshold used to split
+                sharp corners from smooth curve-fit segments.
             high_order_elements: Enable high-order geometric mesh elements.
             high_order_order: Polynomial order for high-order elements.
             high_order_optimize: Run gmsh high-order optimization after meshing.
@@ -1179,6 +1194,7 @@ def mesh(
             curve_fit_layers=curve_fit_layers,
             curve_fit_tolerance_um=curve_fit_tolerance_um,
             curve_fit_min_points=curve_fit_min_points,
+            curve_fit_corner_angle_deg=curve_fit_corner_angle_deg,
             high_order_elements=high_order_elements,
             high_order_order=high_order_order,
             high_order_optimize=high_order_optimize,

src/gsim/palace/base.py

@@ -6,7 +6,7 @@
 import math
 from dataclasses import dataclass, field
 from pathlib import Path
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Literal
 
 import gmsh
 
@@ -220,10 +220,11 @@ def generate_mesh(
     pec_blocks: list[PECBlockConfig] | None = None,
     absorbing_boundary: bool = True,
     merge_via_distance: float = 2.0,
-    curve_fit_mode: str = "line",
+    curve_fit_mode: Literal["line", "spline", "bspline"] = "line",
     curve_fit_layers: list[str] | None = None,
     curve_fit_tolerance_um: float = 0.0,
     curve_fit_min_points: int = 8,
+    curve_fit_corner_angle_deg: float = 45.0,
     high_order_elements: bool = False,
     high_order_order: int = 2,
     high_order_optimize: bool = True,
@@ -256,6 +257,8 @@ def generate_mesh(
         curve_fit_layers: Layer names where curve fitting is applied
         curve_fit_tolerance_um: Point merge tolerance before curve fitting
         curve_fit_min_points: Min contour points required for curve fitting
+        curve_fit_corner_angle_deg: Turn-angle threshold used to identify
+            sharp corners during curve fitting segmentation
         high_order_elements: Enable high-order geometric mesh elements
         high_order_order: Polynomial order for high-order elements
         high_order_optimize: Run gmsh high-order optimization after meshing
@@ -324,6 +327,7 @@ def generate_mesh(
             curve_fit_layers=curve_fit_layers,
             curve_fit_tolerance_um=curve_fit_tolerance_um,
             curve_fit_min_points=curve_fit_min_points,
+            curve_fit_corner_angle_deg=curve_fit_corner_angle_deg,
         )
 
         all_dielectric_tags = {

src/gsim/palace/mesh/generator.py

@@ -542,6 +542,7 @@ def add_patterned_dielectrics(
     curve_fit_layers: list[str] | None = None,
     curve_fit_tolerance_um: float = 0.0,
     curve_fit_min_points: int = 8,
+    curve_fit_corner_angle_deg: float = 45.0,
 ) -> dict[str, list[int]]:
     """Add patterned dielectric volumes from stack dielectric layers.
 
@@ -559,6 +560,8 @@ def add_patterned_dielectrics(
         curve_fit_layers: Layer names where spline/bspline fitting is allowed.
         curve_fit_tolerance_um: Point merge tolerance before curve fitting.
         curve_fit_min_points: Minimum contour points to attempt curve fitting.
+        curve_fit_corner_angle_deg: Turn-angle threshold for corner detection
+            during spline/bspline segmentation.
 
     Returns:
         Dict mapping dielectric layer name -> list of volume tags.
@@ -603,6 +606,7 @@ def add_patterned_dielectrics(
                 loop_mode=surface_loop_mode,
                 fit_tolerance_um=curve_fit_tolerance_um,
                 min_points_for_curve_fit=curve_fit_min_points,
+                corner_turn_threshold_deg=curve_fit_corner_angle_deg,
             )
             if surfacetag is not None:
                 surfaces.append(surfacetag)

src/gsim/palace/mesh/geometry.py

@@ -4,12 +4,15 @@
 
 import logging
 import math
+from itertools import pairwise
 
 import gmsh
 import numpy as np
 
 logger = logging.getLogger(__name__)
 
+_CORNER_TURN_THRESHOLD_DEG = 45.0
+
 # ---------------------------------------------------------------------------
 # Minimalistic meshwell-style entity + boolean pipeline
 # ---------------------------------------------------------------------------
@@ -304,6 +307,7 @@ def _create_wire_loop(
     meshseed: float = 0,
     loop_mode: str = "line",
     point_merge_tol: float = 0.0,
+    corner_turn_threshold_deg: float = _CORNER_TURN_THRESHOLD_DEG,
 ) -> int | None:
     """Create a closed curve loop from polygon vertices.
 
@@ -322,25 +326,105 @@ def _create_wire_loop(
         if math.hypot(x - px, y - py) > tol:
             cleaned.append((x, y))
 
-    if len(cleaned) >= 2:
-        if math.hypot(cleaned[0][0] - cleaned[-1][0], cleaned[0][1] - cleaned[-1][1]) <= tol:
-            cleaned.pop()
+    if len(cleaned) >= 2 and (
+        math.hypot(cleaned[0][0] - cleaned[-1][0], cleaned[0][1] - cleaned[-1][1])
+        <= tol
+    ):
+        cleaned.pop()
 
     if len(cleaned) < 3:
         return None
 
     verts = [kernel.addPoint(x, y, z, meshseed, -1) for x, y in cleaned]
 
+    turn_threshold_deg = max(0.0, min(float(corner_turn_threshold_deg), 180.0))
+
+    def _turn_angle_deg(
+        prev_pt: tuple[float, float],
+        curr_pt: tuple[float, float],
+        next_pt: tuple[float, float],
+    ) -> float | None:
+        ax = curr_pt[0] - prev_pt[0]
+        ay = curr_pt[1] - prev_pt[1]
+        bx = next_pt[0] - curr_pt[0]
+        by = next_pt[1] - curr_pt[1]
+        norm_a = math.hypot(ax, ay)
+        norm_b = math.hypot(bx, by)
+        if norm_a <= tol or norm_b <= tol:
+            return None
+        cos_theta = (ax * bx + ay * by) / (norm_a * norm_b)
+        cos_theta = max(-1.0, min(1.0, cos_theta))
+        return math.degrees(math.acos(cos_theta))
+
+    def _corner_indices(points: list[tuple[float, float]]) -> list[int]:
+        if len(points) < 3:
+            return []
+        corners: list[int] = []
+        for idx in range(len(points)):
+            prev_pt = points[(idx - 1) % len(points)]
+            curr_pt = points[idx]
+            next_pt = points[(idx + 1) % len(points)]
+            turn_angle = _turn_angle_deg(prev_pt, curr_pt, next_pt)
+            if turn_angle is not None and turn_angle >= turn_threshold_deg:
+                corners.append(idx)
+        return corners
+
+    def _corner_to_corner_segments(
+        num_points: int,
+        corners: list[int],
+    ) -> list[list[int]]:
+        unique = sorted(set(corners))
+        if num_points < 3 or len(unique) < 2:
+            return []
+
+        segments: list[list[int]] = []
+        for idx, start in enumerate(unique):
+            end = unique[(idx + 1) % len(unique)]
+            segment = [start]
+            cursor = (start + 1) % num_points
+            while cursor != end:
+                segment.append(cursor)
+                cursor = (cursor + 1) % num_points
+            segment.append(end)
+            segments.append(segment)
+        return segments
+
     if loop_mode in {"spline", "bspline"} and len(verts) >= 4:
+        corners = _corner_indices(cleaned)
+        if len(corners) >= 2:
+            try:
+                segmented_curves: list[int] = []
+                for segment in _corner_to_corner_segments(len(verts), corners):
+                    segment_pts = [verts[i] for i in segment]
+                    if len(segment_pts) >= 4:
+                        if loop_mode == "bspline":
+                            curve = kernel.addBSpline(segment_pts, -1)
+                        else:
+                            curve = kernel.addSpline(segment_pts, -1)
+                        segmented_curves.append(curve)
+                    else:
+                        for p1, p2 in pairwise(segment_pts):
+                            segmented_curves.append(kernel.addLine(p1, p2, -1))
+
+                if len(segmented_curves) >= 3:
+                    return kernel.addCurveLoop(segmented_curves, tag=-1)
+            except Exception:
+                logger.debug(
+                    "Failed segmented %s loop, falling back to full-curve loop",
+                    loop_mode,
+                )
+
         try:
-            curve_pts = verts + [verts[0]]
+            curve_pts = [*verts, verts[0]]
             if loop_mode == "bspline":
                 curve = kernel.addBSpline(curve_pts, -1)
             else:
                 curve = kernel.addSpline(curve_pts, -1)
             return kernel.addCurveLoop([curve], tag=-1)
         except Exception:
-            logger.debug("Failed to create %s loop, falling back to line loop", loop_mode)
+            logger.debug(
+                "Failed to create %s loop, falling back to line loop", loop_mode
+            )
 
     lines = []
     for v in range(len(verts)):
@@ -364,6 +448,7 @@ def create_polygon_surface(
     loop_mode: str = "line",
     fit_tolerance_um: float = 0.0,
     min_points_for_curve_fit: int = 8,
+    corner_turn_threshold_deg: float = _CORNER_TURN_THRESHOLD_DEG,
 ) -> int | None:
     """Create a planar surface from polygon vertices at z height.
 
@@ -382,6 +467,8 @@ def create_polygon_surface(
         fit_tolerance_um: Point merge tolerance before loop creation
         min_points_for_curve_fit: Minimum contour points required to attempt
             spline/bspline loop creation
+        corner_turn_threshold_deg: Turn-angle threshold used to detect corners
+            when segmenting spline/bspline loops
 
     Returns:
         Surface tag, or None if polygon is invalid
@@ -401,6 +488,7 @@ def create_polygon_surface(
         meshseed,
         loop_mode=effective_mode,
         point_merge_tol=fit_tolerance_um,
+        corner_turn_threshold_deg=corner_turn_threshold_deg,
     )
     if exterior_loop is None and effective_mode != "line":
         exterior_loop = _create_wire_loop(
@@ -411,6 +499,7 @@ def create_polygon_surface(
             meshseed,
             loop_mode="line",
             point_merge_tol=fit_tolerance_um,
+            corner_turn_threshold_deg=corner_turn_threshold_deg,
         )
     if exterior_loop is None:
         return None
@@ -432,6 +521,7 @@ def create_polygon_surface(
             meshseed,
             loop_mode="line",
             point_merge_tol=fit_tolerance_um,
+            corner_turn_threshold_deg=corner_turn_threshold_deg,
         )
         if hloop is not None:
             hsurf = kernel.addPlaneSurface([hloop], tag=-1)

src/gsim/palace/mesh/gmsh_utils.py

@@ -53,6 +53,7 @@ class MeshConfig(BaseModel):
     curve_fit_layers: list[str] = Field(default_factory=lambda: ["core", "core2"])
     curve_fit_tolerance_um: float = Field(default=0.0, ge=0)
     curve_fit_min_points: int = Field(default=8, ge=3)
+    curve_fit_corner_angle_deg: float = Field(default=45.0, gt=0, lt=180)
     high_order_elements: bool = False
     high_order_order: int = Field(default=2, ge=2, le=6)
     high_order_optimize: bool = True