r"""Exploration: FBPINN window functions on GMSH-based, non-Cartesian subdomains. The Cartesian FBPINN examples (:mod:`laplacian_1d_fbpinn`, :mod:`laplacian_2d_fbpinn`) decompose a square into a tensor grid of subdomains and write the window functions by hand, in terms of that grid's ``(ix, iy)`` indices. This does not generalize to an arbitrary domain shape. The plan to lift that restriction has four steps: 1. the user creates a domain as usual (a :class:`VolumetricDomain`, e.g. :class:`Square2D`, :class:`Disk2D`, :class:`Flower2D`, ...); 2. GMSH meshes that domain into a small number ``N`` of quadrilateral **macro-cells** -- these become the FBPINN subdomains; 3. a window function is built directly from that :class:`MacroMesh`, one raised-cosine bump per macro-cell, generalizing the tensor-product construction used on Cartesian grids; 4. (not here yet) an FBPINN is trained with these window functions. This file covers steps 1-3, and plots: - the GMSH subdomain partition, for several domain shapes and several (small) numbers of subdomains; - the per-subdomain window functions and their sum, for a handful of subdomains at a time, on a few representative shapes. Step 2 reuses :func:`macro_mesh_from_points` exactly as ``laplacian_2d_fish.py`` already does for FEM (mesh the domain's own boundary curves, obtained from :meth:`VolumetricDomain.full_bc_domain`) for a SIMPLY CONNECTED domain (an open chain of curves, e.g. :class:`Square2D`, :class:`Polygon2D`); :func:`macro_mesh_ogrid_disk` for the one shape (the disk) where an exact, directly-controllable subdomain count is available without GMSH's free mesher; and, for every other shape (one or more independently-CLOSED curves, e.g. :class:`Disk2D`, :class:`Flower2D`, :class:`Annulus2D`), :func:`macro_mesh_from_loops` -- curve by curve, an exact primitive (:class:`ExactCircle`, GMSH's native ``addCircle``) where GMSH has one, a spline fit otherwise. That last case started out as its own GMSH build duplicated here, near-identically to :func:`macro_mesh_from_points`'s own body (same ``gmsh.initialize``/``generate``/``setOrder``/``optimize``/ ``finalize`` sequence, just built from curve loops instead of point-cloud ones) -- moved into ``scimba_jax.mapping.macro_mesh`` itself once the duplication was pointed out, both because that orchestration already existed once and shouldn't exist twice, and because the exact-circle capability is useful to any caller meshing a circular boundary, not just this file. What stayed here is only the one domain-aware step :func:`macro_mesh_from_loops` can't do without importing ``scimba_jax.domains`` (which it deliberately doesn't): recognizing that a given curve *is* a full circle in the first place (:func:`_curve_to_loop_spec`). See :func:`macro_mesh_from_loops`'s own docstring for why the exact-primitive distinction turned out to matter, not just be nicer (a real, measured gap that plain point-cloud spline fitting, even at 800 samples, did not close). The mesh exists ONLY to build the subdomains and their window functions -- step 4's FBPINN sampler draws collocation points from the domain directly (:class:`DomainSampler`, exactly as the Cartesian examples already do), not from the mesh's own footprint. So the actual thing that must hold is the window functions summing to 1 everywhere in the TRUE domain (by its own SDF), not just inside whatever the mesh happens to cover -- which is also why every plot below is masked with ``domain.is_outside``, the domain's own predicate, rather than anything mesh-derived. Step 3's window functions themselves -- :func:`make_mesh_window_function`, one raised-cosine bump per macro-cell, and its Cartesian-grid sibling :func:`make_cartesian_window_function` -- live in the library (``scimba_jax.nonlinear_approximation.approximation_spaces. finite_basis_approximation_spaces``, alongside :class:`FiniteBasisApproximationSpace` itself), since step 4 (an actual FBPINN, see ``laplacian_2d_fbpinn_disk.py``) needs them too and they are not specific to this exploration. See that module's own docstring for the construction itself and the full history of what was tried and measured to be worse (a corner-only straight-line fit left 28% of the O-grid disk's interior with all 5 raw windows simultaneously 0; a per-point Newton-based local frame was smooth in VALUE but had Laplacian spikes past -1500 at cell corners, silently poisoning a Poisson residual). A *coarse* mesh whose cells don't quite reach the domain's true (SDF) boundary is not just a cosmetic gap: right where that gap is widest, every raw window can drop to exactly 0 at once (same 0/0 as above, from discretization error rather than a wrong construction) -- since the window functions must sum to 1 over the TRUE domain (see above), not just inside the mesh, this is a real correctness bug, not a plotting artifact. Measured on the annulus's coarsest kept tier (18 cells, back when it was meshed via a spline through sampled points): a 4.5% mesh-vs-domain area mismatch (``mesh_covers_domain``) and 42 of 8856 interior grid points landing exactly in the resulting undefined band, all within 0.04 of the true boundary (by SDF value). Neither more boundary samples (200 -> 800) nor a higher geometric ``order`` (2 -> 5) moved either number -- the spline was already a visually exact fit to the circle on its own; the gap only appeared once GMSH's OCC kernel built a plane surface bounded by the annulus's two independently-fit splines and meshed it. Switching that closed-curve case to :class:`ExactCircle`'s ``addCircle`` (the actual fix, not a smaller/refined version of the spline) took the SAME-sized tier down to a 0.32% area mismatch and, checked over a finer 150x150 grid, 0 undefined points -- confirmed the same way on every ANNULUS tier kept in this file, all exactly 0. :class:`Flower2D` has no exact GMSH primitive to fall back on (a general polar curve): its coarsest tier (the one actually fed to :func:`make_mesh_window_function` for the window-function plot) also reaches 0 undefined points, but its two finer tiers (kept only to show the partition count growing, never plotted as windows here) still have a few dozen out of ~9200 -- so the annulus's fix does not generalize to "every closed curve is now fine at every tier", only to the one class of curve (an exact circle) GMSH can represent without a fit at all. See the "flower" entry in the ``SHAPES`` catalog below for the (real, measured, not just theoretical) limits of raising the geometric order instead. """ # %% import timeit from pathlib import Path import jax import jax.numpy as jnp import matplotlib.pyplot as plt import numpy as np from scimba_jax.domains.meshless_domains.base import VolumetricDomain from scimba_jax.domains.meshless_domains.domains_2d import Disk2D, Polygon2D, Square2D from scimba_jax.domains.meshless_domains.examples_domains import Annulus2D, Flower2D from scimba_jax.mapping.macro_mesh import ( ExactCircle, MacroMesh, macro_mesh_from_loops, macro_mesh_from_points, macro_mesh_ogrid_disk, ) from scimba_jax.nonlinear_approximation.approximation_spaces.finite_basis_approximation_spaces import ( make_mesh_window_function, ) _HERE = Path(__file__).resolve().parent _OUT = _HERE / "gmsh_subdomains_exploration_output" _OUT.mkdir(exist_ok=True) OVERLAP = 1.5 # same meaning as in laplacian_2d_fbpinn.py: >1 means overlap. N_BOUNDARY_SAMPLES = 200 # %% # ───────────────────────────────────────────────────────────────────────── # Step 1+2: domain -> ordered boundary point cloud -> GMSH MacroMesh. # ───────────────────────────────────────────────────────────────────────── def _sample_curve(curve, n: int, endpoint: bool) -> np.ndarray: """Sample one boundary curve over its own parametric range. Args: curve: A ``SurfacicDomain`` from ``full_bc_domain()``. n: Number of samples. endpoint: Whether to include the curve's own endpoint (only makes sense for a curve that is a full loop by itself). Returns: Physical points ``(n, 2)``. """ t0, t1 = (float(v) for v in curve.parametric_domain.bounds[0]) ts = jnp.linspace(t0, t1, n, endpoint=endpoint)[:, None] return np.asarray(curve.surface_o_mapping(ts)) def _is_closed_curve(curve, atol: float = 1e-6) -> bool: """Whether a single curve already closes on itself (e.g. a full circle). Args: curve: A ``SurfacicDomain`` from ``full_bc_domain()``. atol: Absolute tolerance on the start/end physical distance. Returns: True if the curve's image at its two parametric endpoints coincide. """ t0, t1 = (float(v) for v in curve.parametric_domain.bounds[0]) p0 = np.asarray(curve.surface_o_mapping(jnp.array([[t0]])))[0] p1 = np.asarray(curve.surface_o_mapping(jnp.array([[t1]])))[0] return float(np.linalg.norm(p1 - p0)) < atol def _signed_area(loop: np.ndarray) -> float: """Shoelace signed area of an ordered closed polygon. Args: loop: Ordered points ``(n, 2)`` (first point not repeated). Returns: The signed area (positive if CCW). """ x, y = loop[:, 0], loop[:, 1] return 0.5 * float(np.sum(x * np.roll(y, -1) - np.roll(x, -1) * y)) def domain_boundary_curves( domain: VolumetricDomain, n_per_curve: int = 400 ) -> list[np.ndarray]: """Every closed loop of a domain's TRUE boundary, for plotting. The window-function plots are masked with ``domain.is_outside`` (the domain's own SDF), so they are already evaluated on the true domain, not the meshed one -- but a masked ``contourf`` on a finite grid still shows a jagged, staircase-like edge that can look mesh-like at a glance. This draws the exact analytic boundary curve(s) on top (one loop for a simply connected domain, several for a multiply connected one like :class:`Annulus2D`), so the true, smooth boundary is unambiguous regardless of grid resolution -- reusing the same closed/open curve classification :func:`mesh_domain` uses to build the mesh in the first place, just for plotting instead of meshing. Args: domain: The volumetric domain. n_per_curve: Samples per boundary curve. Returns: One or more ordered, closed point loops (first point not repeated). """ curves = domain.full_bc_domain() if any(_is_closed_curve(c) for c in curves): return [_sample_curve(c, n_per_curve, endpoint=False) for c in curves] return [ np.concatenate( [_sample_curve(c, n_per_curve, endpoint=False) for c in curves], axis=0 ) ] def domain_boundary_loop( domain: VolumetricDomain, n_per_curve: int = N_BOUNDARY_SAMPLES ) -> np.ndarray: """Ordered, closed boundary point cloud for a SIMPLY CONNECTED domain. Generalizes ``fish_boundary_points`` from ``laplacian_2d_fish.py`` to any :class:`VolumetricDomain` whose boundary is several open curves that chain head-to-tail into one loop (:class:`Square2D`'s 4 edges, :class:`Polygon2D`'s edges, ...): each curve is sampled over its own parametric range with ``endpoint=False``, so consecutive curves join without duplicating the shared corner and the whole thing closes automatically -- exactly what :func:`macro_mesh_from_points` expects. A domain whose boundary is already one or more independently-closed curves (:class:`Disk2D`, :class:`Flower2D`, :class:`Annulus2D`, ...) does NOT go through this path -- see :func:`mesh_domain`. Args: domain: The volumetric domain to mesh. n_per_curve: Number of samples per boundary curve. Returns: The boundary point cloud, in the format :func:`macro_mesh_from_points` expects. Raises: ValueError: If any curve of ``full_bc_domain()`` is already closed on its own -- that is :func:`mesh_domain`'s other path. """ curves = domain.full_bc_domain() if any(_is_closed_curve(c) for c in curves): raise ValueError( "domain_boundary_loop() is for an open chain of curves; " "a closed curve should go through mesh_domain()'s other path." ) pts = [_sample_curve(c, n_per_curve, endpoint=False) for c in curves] return np.concatenate(pts, axis=0) def _curve_to_loop_spec(curve, sampled: np.ndarray) -> ExactCircle | np.ndarray: """Classify one closed boundary curve for :func:`macro_mesh_from_loops`. The only domain-aware step left in this file for the closed-curve case: is this curve a full :class:`Circle2D` (:class:`Disk2D`'s own, and each of :class:`Annulus2D`'s two), which GMSH can represent exactly? All the GMSH orchestration itself -- including the ``addCircle`` vs spline-fit choice this feeds into -- lives in :func:`macro_mesh_from_loops`, not here; see its docstring for why the distinction turned out to matter, not just be nicer (a real, measured mesh-vs-domain area gap that neither more samples nor a higher geometric order closed). Args: curve: A closed ``SurfacicDomain`` (a full loop by itself). sampled: This curve's own sampled point cloud, already needed by the caller to pick the outer loop by enclosed area -- reused here instead of re-sampling, and simply ignored for a full circle. Returns: An :class:`ExactCircle` for a full circle, else ``sampled`` as-is. """ if curve.domain_type == "Circle2D": return ExactCircle( float(curve.center[0]), float(curve.center[1]), float(curve.radius) ) return sampled def mesh_domain( domain: VolumetricDomain, mesh_size: float, order: int = 2, smooth: bool = True, n_per_curve: int = N_BOUNDARY_SAMPLES, ) -> MacroMesh: """Step 2: mesh a domain's own boundary into a small-``N`` :class:`MacroMesh`. Two paths, chosen from ``domain.full_bc_domain()``'s own curves: - simply connected (an open chain of curves, e.g. :class:`Square2D`, :class:`Polygon2D`): the point-cloud + :func:`macro_mesh_from_points` path (already exact for straight edges, via ``smooth=False``). - one or more independently-closed curves (:class:`Disk2D`, :class:`Flower2D`, :class:`Annulus2D`, ...): :func:`macro_mesh_from_loops`, curve by curve, using an exact primitive where GMSH has one (a full circle, via :func:`_curve_to_loop_spec`) and a spline fit otherwise. Args: domain: The volumetric domain to mesh. mesh_size: GMSH target characteristic length (coarser -> fewer, bigger subdomains). order: Polynomial degree of the element geometry. smooth: Fit a spline through the boundary (curved shapes) or join with straight segments (polygonal shapes, to keep sharp corners exact). Ignored for a closed curve with an exact primitive (always exact then, regardless of ``smooth``). n_per_curve: Boundary samples per curve fed to GMSH (fallback spline case only). Returns: The assembled :class:`MacroMesh`, one macro-cell per subdomain. """ curves = domain.full_bc_domain() if not any(_is_closed_curve(c) for c in curves): # A straight edge only needs its own endpoint: with smooth=False, # macro_mesh_from_points() joins consecutive points with straight # lines, so sampling 200 points per edge would turn each polygon # side into 200 tiny collinear segments instead of one straight line. effective_n_per_curve = n_per_curve if smooth else 1 boundary = domain_boundary_loop(domain, effective_n_per_curve) return macro_mesh_from_points( boundary, order=order, mesh_size=mesh_size, smooth=smooth ) sampled = [_sample_curve(c, n_per_curve, endpoint=False) for c in curves] outer_idx = ( int(np.argmax([abs(_signed_area(loop)) for loop in sampled])) if len(curves) > 1 else 0 ) specs = [_curve_to_loop_spec(c, loop) for c, loop in zip(curves, sampled)] holes = [s for k, s in enumerate(specs) if k != outer_idx] return macro_mesh_from_loops( specs[outer_idx], holes=holes or None, order=order, mesh_size=mesh_size, smooth=smooth, ) # %% # ───────────────────────────────────────────────────────────────────────── # Plotting. # ───────────────────────────────────────────────────────────────────────── def _cell_outline(mesh: MacroMesh, cell_idx: int, n_edge: int = 8) -> np.ndarray: """Physical outline of one macro-cell, sampling its (possibly curved) edges. Args: mesh: The macro-mesh. cell_idx: The macro-cell index. n_edge: Samples per edge. Returns: Ordered physical boundary points of the cell, ``(4 * n_edge, 2)``. """ t = np.linspace(0.0, 1.0, n_edge, endpoint=False) ref = np.concatenate( [ np.stack([t, np.zeros_like(t)], axis=1), np.stack([np.ones_like(t), t], axis=1), np.stack([1.0 - t, np.ones_like(t)], axis=1), np.stack([np.zeros_like(t), 1.0 - t], axis=1), ] ) phys = jax.vmap(lambda p: mesh.forward(cell_idx, p))(jnp.asarray(ref)) return np.asarray(phys) def mesh_covers_domain( domain: VolumetricDomain, mesh: MacroMesh, n_monte_carlo: int = 200_000, seed: int = 0, is_inside_fn=None, ) -> tuple[float, float, float]: """Sanity check: does the mesh's area match the domain's true area? Catches a real failure mode of :func:`macro_mesh_from_points`: when ``mesh_size`` is coarser than a hole's own feature size, GMSH's mesher can silently mesh straight across the hole instead of excluding it (a node ends up sitting exactly at the disk center) -- the resulting :class:`MacroMesh` is geometrically valid (conforming, no gaps) but wrong, and nothing about it raises. Comparing the summed macro-cell area to a Monte-Carlo estimate of the domain's own area (from ``domain.is_inside``, so this needs nothing beyond what's already used for step 1) catches that generically, for any shape. Args: domain: The domain that was meshed. mesh: The resulting macro-mesh. n_monte_carlo: Number of random samples for the domain-area estimate. seed: RNG seed for the sample points. is_inside_fn: A pre-jitted ``domain.is_inside`` to reuse across several calls for the same domain (a fresh ``jax.jit`` wrapper recompiles from scratch every time, and this is typically called once per mesh_size on the SAME domain). Defaults to jitting ``domain.is_inside`` locally, for standalone use. Returns: ``(mesh_area, domain_area, relative_error)``. """ if is_inside_fn is None: is_inside_fn = jax.jit(domain.is_inside) rng = np.random.default_rng(seed) b = np.asarray(domain.bounds) lo, hi = b[:, 0], b[:, 1] pts = rng.uniform(lo, hi, size=(n_monte_carlo, 2)) # jit, not just vmap: is_inside is already internally vmapped # (DomainMapping batches every pointwise op it wraps), but vmap alone # still dispatches each op in the trace eagerly, one XLA call per op; # jit fuses the whole batched graph into a single compiled call, which # matters a lot at 200k points. inside = np.asarray(is_inside_fn(jnp.asarray(pts))).reshape(-1) domain_area = float(np.prod(hi - lo)) * float(inside.mean()) mesh_area = sum( abs(_signed_area(_cell_outline(mesh, i, n_edge=6))) for i in range(mesh.n_cells) ) rel_err = abs(mesh_area - domain_area) / domain_area return mesh_area, domain_area, rel_err def plot_partition(mesh: MacroMesh, ax, title: str) -> None: """Plot one macro-mesh's subdomain partition on a given axis. Args: mesh: The macro-mesh to draw. ax: A matplotlib axis. title: Axis title (the subdomain count is appended). """ cmap = plt.get_cmap("tab20") for i in range(mesh.n_cells): outline = _cell_outline(mesh, i) ax.fill( outline[:, 0], outline[:, 1], color=cmap(i % 20), alpha=0.6, edgecolor="k", linewidth=0.6, ) ax.set_aspect("equal") ax.set_title(f"{title}\n({mesh.n_cells} subdomains)") def plot_window_functions( domain: VolumetricDomain, mesh: MacroMesh, overlap: float, title: str, n_grid: int = 180, ): """Plot each subdomain's window function and their sum (step 3, "low N"). Mirrors the ``plot_window_functions`` block of ``laplacian_2d_fbpinn.py``, but driven by the mesh's macro-cells instead of a hardcoded Cartesian grid, and evaluated on the domain's TRUE (SDF) interior, not the meshed one -- the mesh only exists to build the subdomains/window functions in the first place (see the module docstring), it does not bound where a future step 4's sampler will draw points. Concretely: every point is masked with ``domain.is_outside`` (masking with the mesh's own footprint would be the wrong thing here, even though the two are close), and the domain's own exact boundary curve(s) are drawn on top of every subplot (:func:`domain_boundary_curves`) so the true, smooth boundary is unambiguous regardless of the (still finite, still slightly jagged at its edge) contour grid. Args: domain: The domain that was meshed (for masking outside the shape, and for drawing its true boundary on top of every subplot). mesh: The macro-mesh whose cells are the subdomains. overlap: Passed to :func:`make_mesh_window_function`. title: Figure title. n_grid: Grid resolution per axis for the contour plots. Returns: The created matplotlib figure. """ _window_function, n_subdomains, all_windows = make_mesh_window_function( mesh, overlap ) b = domain.get_extended_bounds(factor=0.05) xs = jnp.linspace(b[0, 0], b[0, 1], n_grid) ys = jnp.linspace(b[1, 0], b[1, 1], n_grid) xx, yy = jnp.meshgrid(xs, ys, indexing="ij") pts = jnp.stack([xx.ravel(), yy.ravel()], axis=-1) outside = np.asarray(jax.jit(domain.is_outside)(pts)).reshape(-1) boundary_loops = domain_boundary_curves(domain) fig, axes = plt.subplots( 1, n_subdomains + 1, figsize=(3.2 * (n_subdomains + 1), 3.4), squeeze=False ) axes = axes[0] # jit the batched all_windows (O(n_subdomains) raw bumps per point), NOT # n_subdomains separate jits of window_function (O(n_subdomains) raw # bumps EACH, i.e. O(n_subdomains^2) total, and n_subdomains recompiles): # measured on the annulus's 18-subdomain mesh, the naive # "jit(window_function) inside the loop" version took 23s to compile for # this reason alone. jit still matters (vs plain vmap): is_inside's note # above applies the same way here. windows_by_point = jax.jit(jax.vmap(all_windows))(pts) # (n_points, n_subdomains) sum_windows = jnp.sum(windows_by_point, axis=1) def _draw_true_boundary(ax) -> None: for loop in boundary_loops: closed = np.concatenate([loop, loop[:1]], axis=0) ax.plot(closed[:, 0], closed[:, 1], color="black", linewidth=0.8) for i in range(n_subdomains): w = windows_by_point[:, i] w_plot = np.where(outside, np.nan, np.asarray(w)) axes[i].contourf(xx, yy, w_plot.reshape(n_grid, n_grid), levels=20) _draw_true_boundary(axes[i]) axes[i].set_title(f"window {i}") axes[i].set_aspect("equal") sw_plot = np.where(outside, np.nan, np.asarray(sum_windows)) im = axes[-1].contourf(xx, yy, sw_plot.reshape(n_grid, n_grid), levels=20) _draw_true_boundary(axes[-1]) plt.colorbar(im, ax=axes[-1]) axes[-1].set_title("sum of windows") axes[-1].set_aspect("equal") fig.suptitle(title) fig.tight_layout() return fig # %% # ───────────────────────────────────────────────────────────────────────── # Several domain shapes, several (small) numbers of subdomains. # ───────────────────────────────────────────────────────────────────────── # Straight-edged (polygonal) shapes are meshed with smooth=False, so GMSH # joins the sampled boundary points with straight segments and keeps corners # exact instead of rounding them off with a spline fit. # # N is only *indirectly* controlled through mesh_size (GMSH's free mesher # decides how many cells actually fit), and does so in discrete, non-monotone # jumps rather than smoothly -- e.g. the square sits at exactly 24 cells for # every mesh_size in [0.5, 0.9]. So mesh_size is picked per shape, by probing # a few values and keeping three that give a clearly increasing, low-to- # moderate subdomain count, rather than one global rule (a fixed fraction of # the bounding-box diagonal) that happened to alias to the same count twice # for the square and the L-shape. SHAPES: dict[str, tuple[VolumetricDomain, bool, list[float], int]] = { "square": ( Square2D([(-1.0, 1.0), (-1.0, 1.0)], is_main_domain=True), False, [1.131, 0.9, 0.3], 2, ), "l_shape": ( Polygon2D( [(-1, -1), (1, -1), (1, 0), (0, 0), (0, 1), (-1, 1)], is_main_domain=True, ), False, [1.131, 0.7, 0.3], 2, ), "disk": ( Disk2D([0.0, 0.0], 1.0, is_main_domain=True), True, [0.884, 0.566, 0.354], 2, ), # No exact GMSH primitive for a general polar curve (unlike the circles # below), so raising the geometric order was tried as the remaining # lever: order=4 measurably shrinks the mesh-vs-domain area error on the # coarsest tier (0.95% -> ~0.2-0.3%) -- but is NOT SAFE here. GMSH's # HighOrderOptimize pass can hard-abort (a C++ exception outside # Python's try/except, killing the whole process) with "Failed to reach # critical value ... ScaledJac", and it did so NON-DETERMINISTICALLY: the # SAME (mesh_size=1.193, order=4) call, in a fresh process, succeeded # some runs and aborted others. Order=3 and order=4 both aborted at least # once across the three tiers kept here. The flower's tight petal # "waists" (a=0.35 relative to R=1.0) are apparently right at the edge of # what a coarse curved boundary can represent with a valid (positive) # Jacobian everywhere, and whatever GMSH's optimizer uses to search for # one isn't reproducible run to run. Kept at order=2 -- slower to # converge, but never once crashed across dozens of runs, and (checked # the way that actually matters -- see the module docstring) already # gives 0 undefined window-function points on the COARSEST tier, the one # the window-function plot actually uses. The other two, finer tiers # (kept only to show the partition count growing, never fed to # make_mesh_window_function() in this file) do still have a few dozen such # points out of ~9200 checked -- consistent with GMSH's free # recombination not improving monotonically with refinement on this # shape (see the "flower" entry in the mesh_size sweep printout: area # error goes UP from the 11- to the 22-cell tier before coming back # down) -- a real caveat for whichever tier a future step 4 picks, not # something to silently round up to "fixed everywhere". "flower": ( Flower2D(R=1.0, a=0.35, n=5, is_main_domain=True), True, [1.193, 0.764, 0.477], 2, ), # Before mesh_domain() used an exact ExactCircle/addCircle for a # full-circle boundary (see macro_mesh_from_loops), a mesh_size coarser # than the hole's own radius # (0.35) didn't just under-resolve the hole -- it made GMSH mesh straight # across it (a node landing exactly at the domain center): mesh_size=0.884 # or 0.566 both silently produced a plain, hole-free disk. Caught by # mesh_covers_domain() below. Gone with the exact circle (checked up to # mesh_size=1.5, comfortably coarser than the domain itself), which is # also why these three are much coarser now than that fix needed. With # an exact primitive, order stops mattering here (checked: order=2 vs 4 # barely moves the already-small error), so it's kept at 2. "annulus": ( Annulus2D(0.35, 1.0, is_main_domain=True), True, [0.9, 0.6, 0.354], 2, ), } print( "\n\n@@@@@@@@@@@@@@@ step 1+2: domain -> GMSH macro-mesh @@@@@@@@@@@@@@@@@@@@@@\n" ) meshes: dict[str, list[MacroMesh]] = {} for name, (domain, smooth, mesh_sizes, order) in SHAPES.items(): meshes[name] = [] is_inside_jit = jax.jit( domain.is_inside ) # one compile, reused across mesh_sizes below fig, axes = plt.subplots(1, len(mesh_sizes), figsize=(4.2 * len(mesh_sizes), 4.2)) for ax, mesh_size in zip(axes, mesh_sizes): start = timeit.default_timer() mesh = mesh_domain(domain, mesh_size, order=order, smooth=smooth) end = timeit.default_timer() meshes[name].append(mesh) mesh_area, domain_area, rel_err = mesh_covers_domain( domain, mesh, is_inside_fn=is_inside_jit ) print( f"{name:>10s}: mesh_size={mesh_size:.4f} -> {mesh.n_cells:3d} subdomains, " f"{mesh.n_interfaces:3d} interfaces, {mesh.n_boundary_faces:3d} boundary faces, " f"area {mesh_area:.4f} vs domain {domain_area:.4f} " f"(rel. err {rel_err:.2%}) ({end - start:.2f} s)" ) if rel_err > 0.08: print( f" ! mesh area is off by {rel_err:.1%}: the mesh likely does not " "cover the domain correctly (e.g. mesh_size coarser than a hole's " "feature size) -- do not trust this partition." ) plot_partition(mesh, ax, f"{name}, h={mesh_size:.3f}") fig.suptitle(f"GMSH subdomain partitions of the '{name}' domain") fig.tight_layout() fig.savefig(_OUT / f"partition_{name}.png", dpi=150) # The disk also gets the exact-N O-grid alternative (block-structured, no # free GMSH meshing needed): n=1 -> 5 cells, n=2 -> 20, n=3 -> 45 (5*n^2). print() disk_domain, _, _, _ = SHAPES["disk"] fig, axes = plt.subplots(1, 3, figsize=(4.2 * 3, 4.2)) for ax, n in zip(axes, (1, 2, 3)): ogrid_mesh = macro_mesh_ogrid_disk(radius=1.0, inner=0.5, n=n) print(f" disk (O-grid, n={n}): {ogrid_mesh.n_cells:3d} subdomains (exact)") plot_partition(ogrid_mesh, ax, f"O-grid disk, n={n}") fig.suptitle("Disk: exact subdomain count via the block-structured O-grid") fig.tight_layout() fig.savefig(_OUT / "partition_disk_ogrid.png", dpi=150) # %% # ───────────────────────────────────────────────────────────────────────── # Step 3: window functions, for the coarsest (smallest-N) mesh per shape. # ───────────────────────────────────────────────────────────────────────── print( "\n\n@@@@@@@@@@@@@@@ step 3: window functions on the coarsest mesh @@@@@@@@@@@@@@@@@@@@@@\n" ) for name, (domain, _smooth, _mesh_sizes, _order) in SHAPES.items(): coarsest = meshes[name][0] print(f"{name:>10s}: plotting {coarsest.n_cells} window functions") fig = plot_window_functions( domain, coarsest, OVERLAP, title=f"FBPINN window functions on '{name}' ({coarsest.n_cells} subdomains)", ) fig.savefig(_OUT / f"windows_{name}.png", dpi=150) # a lower-N view of the O-grid disk too (n=1 -> 5 cells): the cleanest case # to check the construction against, since the 5 macro-cells are the classic # center-square + 4-petals O-grid decomposition. ogrid_5 = macro_mesh_ogrid_disk(radius=1.0, inner=0.5, n=1) fig = plot_window_functions( disk_domain, ogrid_5, OVERLAP, title="FBPINN window functions on the O-grid disk (5 subdomains)", ) fig.savefig(_OUT / "windows_disk_ogrid.png", dpi=150) print(f"\nAll figures saved under {_OUT}") plt.show() # %%