diff --git a/.claude/skills/adapt-on-top-faults/SKILL.md b/.claude/skills/adapt-on-top-faults/SKILL.md new file mode 100644 index 00000000..a209d5b1 --- /dev/null +++ b/.claude/skills/adapt-on-top-faults/SKILL.md @@ -0,0 +1,236 @@ +--- +name: adapt-on-top-faults +description: Recipe for Underworld3 FAULT models on an NVB adapt-on-top mesh — resolve a fault Surface by LOCAL refinement (mesh.adapt(engine="nvb") returns a child), drive it from the fault's EXACT signed distance, use rotated strong free-slip (composes with transverse-isotropy where Nitsche does not), recover dynamic topography from the constraint reaction, and run advection-diffusion on the adapted mesh with field transfer across re-adaptation. Reach for THIS for instantaneous/coupled fault-flow problems. For MMPDE node-movement convection use the `adaptive-meshing` skill instead; for rendering use `uw-visualisation`. +--- + +# adapt-on-top-faults + +The validated recipe for **fault problems on a locally-refined (adapt-on-top) mesh**. +Distilled from the annulus fault study (2026-07, `feature/adapt-on-top`). + +**This is the REFINEMENT paradigm**, not the mover one: +- `mesh.adapt(metric, engine="nvb")` bisects the base finest **locally** and returns + a **new child mesh** (`child.parent is mesh`). It is *adapt / re-adapt*, NOT node + movement — non-cumulative (each call re-marks from the static base). The child owns + a custom-P geometric-MG (FMG) tail so solvers on it get multigrid for free. +- For **MMPDE / equidistribution node movement** (deforming the same mesh to a field) + use the **`adaptive-meshing`** skill instead. Different tool; don't mix them up. + +Reference implementations (copy from these — all run np1/2/4): +`~/+Simulations/nvb_parallel_fault_study/` (weak-fault Stokes + FMG), +`~/+Simulations/shear_box_fault_study/` (iso vs TI, orientation sweeps), +`~/+Simulations/annulus_fault_study/` (rotated free-slip + topography + moving fault + +advection-diffusion). Companion skills: `uw-visualisation`, `adaptive-meshing`. + +--- + +## The core loop (fault → metric → adapt → child) + +```python +import underworld3 as uw, numpy as np, sympy + +# Base mesh MUST be built with refinement>=1 (supplies the coarse MG tail that NVB +# extends). Base cellSize ~ 2x the target h_near is the SWEET SPOT (see gotchas). +base = uw.meshing.Annulus(radiusInner=0.55, radiusOuter=1.0, cellSize=0.08, + refinement=1, qdegree=3) # or UnstructuredSimplexBox(..., refinement=2) + +# fault as a Surface (polyline control points, N x 3 with z=0 in 2D) +fault = uw.meshing.Surface("fault", base, fault_pts, symbol="F") +fault.discretize() + +# METRIC = a CALLABLE built from the fault's EXACT signed distance. This is the key +# to clean, non-patchy grading: it is evaluated at each refined level's centroids, +# so it resolves itself at the new resolution (no P1-field aliasing). +metric = fault.refinement_metric_function(h_near=0.02, h_far=0.08, width=0.05, + profile="linear") +child = base.adapt(metric, max_levels=3, engine="nvb") # -> graded child mesh +``` + +- `engine="nvb"` = graded newest-vertex bisection (bounded closure, parallel via the + native `uwnvb` transform; bit-confluent serial↔parallel). `engine="sbr"` = uniform + patch (the default; not graded). NVB is 2D only for now. +- `max_levels` is the isotropic-equivalent depth (NVB runs `2*max_levels` bisection + passes). The metric shape decides the grading; `max_levels` just caps it. + +### Metric options (all accepted by `adapt`) +1. **callable** `metric(centroids)->M` — **preferred for faults**. Evaluated per level. +2. MeshVariable / sympy expression — sampled via `uw.function.evaluate` from the BASE + mesh → a peaked `M=1/h²` aliases → *patchy* levels. Avoid for thin features. + +**Custom refinement shape** — pass any callable. For a *fat, uniformly-fine* band +(not just a thin line at the fault), a flat-core metric: +```python +def metric(pts, _f=fault): + d = _f.unsigned_distance(pts) # EXACT distance at arbitrary points + core, ramp, hn, hf = 0.02, 0.05, 0.01, 0.08 + h = np.where(d < core, hn, np.minimum(hn + (hf-hn)*(d-core)/ramp, hf)) + return 1.0 / h**2 +``` + +`Surface` distance API (all exact, arbitrary query points): +`fault.signed_distance(coords)`, `fault.unsigned_distance(coords)`, +`fault.director` (unit normal = normalised ∇(signed distance); the TI weak-plane +director), `fault.refinement_metric_function(...)`. + +--- + +## Fault as a constitutive weak zone (iso and TI) + +The metric only needs the fault GEOMETRY (pure distance). The constitutive weak zone +needs the fault's distance/normal ON THE CHILD. Two ways: + +```python +# (A) re-home the SAME fault onto the child (cleanest; distance recomputes on child) +fault.remap_to(child) +eta = fault.influence_function(width=0.04, value_near=1e-3, value_far=1.0, + profile="smoothstep") # isotropic weak zone + +# (B) or build a child-side Surface (needed if the base fault is still in use for a +# base-mesh metric — symbol disambiguation refuses a base-mesh symbol in a child +# solver). fault_c = uw.meshing.Surface("fault_c", child, fault_pts); fault_c.discretize() +``` + +Isotropic weak zone: +```python +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = eta # drops near fault +``` + +Transverse-isotropic weak PLANE (the physical fault; low fault-parallel shear): +```python +stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_bulk # normal viscosity (constant) +stokes.constitutive_model.Parameters.shear_viscosity_1 = eta # weak near fault -> bulk far +stokes.constitutive_model.Parameters.director = fault.director # unit fault normal +``` +The **TI Jacobian is the full consistent tangent** (not isotropic + defect +correction — that framing is WRONG). The TI velocity FMG V-cycle needs a few more +iters than iso (~8 vs ~2) because of a directional near-null mode the isotropic +point-smoother doesn't damp — bounded and contrast-independent, not a bug. TI needs +~2 elements across the weak zone to resolve (iso ~1). + +--- + +## Rotated strong free-slip (the reason to use this, not Nitsche) + +Nitsche free-slip is INCOMPATIBLE with the TI model (its penalty scales by an +isotropic viscosity). Rotated strong free-slip imposes `u·n̂=0` as an ESSENTIAL +constraint in a per-node (n,t) frame → machine-zero leakage AND composes with TI. + +```python +nhat = mesh.CoordinateSystem.unit_e_0 # exact radial normal (annulus/sphere) +stokes.add_rotated_freeslip_bc("Upper", normal=nhat) +stokes.add_rotated_freeslip_bc("Lower", normal=nhat) +stokes.petsc_use_pressure_nullspace = True # enclosed -> pressure gauge +stokes.solve() # rigid-rotation gauge auto-removed +# Convergence status: read stokes._rotated_freeslip_info = {ksp_reason, +# nonlinear_iterations, rotation_gauge_removed, reaction} — NOT s.snes.getConvergedReason() +# (the rotated solve is a manual loop, not snes.solve). Also sanity-check v·n leakage (~1e-16). +``` + +**Nonlinear rheology / warm-start / timestepping (PR #298, `feature/rotated-snes`):** +rotated free-slip now works *inside* the nonlinear iteration — a nonlinearity probe +auto-dispatches power-law / VEP / TI-with-yield models to a Newton/Picard loop +(`solve_rotated_freeslip_nonlinear`), so warm-started time loops are correct. On the +worktree state *before* #298 lands, the rotated path is a SINGLE linear solve — it +silently returns one Newton linearisation from `u=0` for a nonlinear model. If you +run nonlinear TI + timestepping with rotated free-slip, make sure #298 is in. + +**Dynamic topography** from the constraint reaction (the reason to bother): +```python +h = uw.discretisation.MeshVariable("h", mesh, 1, degree=1) +stokes.dynamic_topography("Upper", h, buoyancy_scale=rho_g) # h = -(σ_nn - mean)/ρg +xs, sig = stokes.boundary_normal_traction("Upper") # or the raw σ_nn (lumped-mass) +``` + +--- + +## FMG under rotated free-slip — READ THIS (rough edge) + +`rotated_bc.solve_rotated_freeslip` builds its OWN fieldsplit KSP and uses custom-P +FMG **only if `solver._custom_mg` is set** (via `set_custom_fmg`). It does **not** +read the native `dm_hierarchy`. So: +- On an `adapt()` **child**, the mesh-owned custom-P tail is auto-picked-up for the + Stokes velocity block → FMG for free (standard, non-rotated solves). +- On a plain **refined `Annulus`** with **rotated** free-slip, native FMG is ignored; + to get FMG you must build a coarse-mesh tail and call: + ```python + from underworld3.utilities.custom_mg import set_custom_fmg + set_custom_fmg(stokes, [Annulus(cs=0.16), Annulus(cs=0.08)], field_id=0) # velocity block + ``` + (~4 velocity iters vs ~26 GAMG). Otherwise it falls back to GAMG (fine, just slower). + +The default velocity-block preconditioner and the pressure Schur (`1/η`) are already +near-optimal for TI — do **not** hand-roll a "TI-aware Schur"; measured, the default +`1/η₀` beats every alternative (the weak TI mode is fault-parallel shear, which is +volume-preserving, so pressure sees η₀). + +--- + +## Moving fault + re-adaptation + field transfer + +`adapt()` is non-cumulative and deterministic. Move the fault, re-adapt from the +static base, carry any field by interpolation: + +```python +for step in range(N): + fault_pts = move(fault_pts, step) # kinematics + fault = uw.meshing.Surface(f"fault{step}", base, fault_pts, symbol="F"); fault.discretize() + child = base.adapt(fault.refinement_metric_function(...), max_levels=3, engine="nvb") + # verify: folded=0 (all cell |vol|>0), base unchanged (non-cumulative) + # carry a field child_{k-1} -> child_k by interpolation: + T = uw.discretisation.MeshVariable(f"T{step}", child, 1, degree=1) + if T_prev is None: + T.data[:,0] = uw.function.evaluate(T0_expr, T.coords) + else: + T.data[:,0] = uw.function.evaluate(T_prev.sym, T.coords) # mesh->mesh interp + T_prev = T +``` +Transfer error is **non-accumulating** for a smooth field (bounded by the per-mesh P1 +representation floor, ~0.2% on a fine base, ~3% on a coarse base) — repeated +re-meshing does not diffuse a smooth field away. Sharp fronts lose more per transfer. + +--- + +## Advection-diffusion on the adapted mesh + +```python +T = uw.discretisation.MeshVariable("T", child, 1, degree=2) +adv = uw.systems.AdvDiffusionSLCN(child, u_Field=T, V_fn=stokes.u.sym, order=1, + monotone_mode="clamp") # clamp = bounded, no overshoot +adv.constitutive_model = uw.constitutive_models.DiffusionModel +adv.constitutive_model.Parameters.diffusivity = 2e-4 +adv.f = 0.0 +dt = 0.015 # FIXED dt — SLCN is semi-Lagrangian (unconditionally + # stable). estimate_dt() reports the tiny fault-band + # Courant limit; do NOT use it to size the step. +for step in range(nsteps): + adv.solve(timestep=dt) +``` +- **SLCN dt is NOT Courant-limited** by the fine fault-band cells — pick dt from the + coarse-region advection, verify accuracy. +- The scalar AD auto-FMG-injection bug (velocity-block custom-P mismatching a scalar + operator on adapt children → PtAP error 60) is **FIXED**: `auto_inject_custom_mg` + now calls `snes.setUp()` before reading the finest reduced map (so the DM section + is the finalized space the operator lives on) and validates that map against the + assembled operator. Custom-P installs successfully on scalar SLCN AdvDiffusion on + NVB adapt children (no skip-guard, no workaround) — `bugfix/custom-mg-parallel`. + +--- + +## Gotchas / rough edges (candidates to fix as we go) + +| symptom / edge | cause & handling | +|---|---| +| patchy along-fault refinement (level 4 here, 2 there) | P1-interpolated `M=1/h²` aliasing. Use the **callable exact-distance** metric. | +| child solver rejects `fault.distance.sym` (foreign-mesh error) | symbol disambiguation. `fault.remap_to(child)` or build a child-side `Surface`. **Rough edge**: needing two Surface objects. | +| adapt added-cell count jitters ±30% step-to-step | small-number geometry on a thin band. Deterministic; quality (near-fault h) is constant. Base ≈ **2× target** minimises it (CoV ~10% at 0.08 base vs 19% at 0.14, 24% at 0.05). | +| tiny AD timestep / slow advection | `estimate_dt` returns the fault-band Courant limit. Use a fixed dt (SLCN is unconditionally stable). **Rough edge**: `estimate_dt` not adapt-aware. | +| surface-breaking fault: huge local vmax; topo peak keeps growing | real stress singularity at the outcrop. Topo peak **saturates** (integrable/log, bounded by finite buoyancy) — far field is fine; only the pointwise outcrop value is mesh-dependent. | +| rotated free-slip Stokes "not converged" | `s.snes` isn't the solving object (manual loop). Read `stokes._rotated_freeslip_info['ksp_reason']` / `['nonlinear_iterations']` (PR #298); sanity-check v·n leakage. | +| nonlinear TI/VEP + rotated free-slip + timestepping gives a wrong (frozen) answer | pre-#298 the rotated path is ONE linear solve (one Newton step from u=0). PR #298 runs it inside a Newton/Picard loop — ensure it's merged for nonlinear/warm-start runs. | +| FMG not used under rotated free-slip | rotated_bc uses a self-contained KSP and reads `solver._custom_mg`, not the native `dm_hierarchy` — `set_custom_fmg(..., field_id=0)`. FUNDAMENTAL (DM-less rotated operator can't use the DM-coupled fieldsplit; #298 keeps it), not a quick fix. | +| NVB at np>1 raises NotImplementedError | native `_nvb_transform` extension not built (needs the custom-PETSc/amr env). | + +**Build/run**: this lives on the `feature/adapt-on-top` worktree; env +`.pixi/envs/amr-dev/bin/{python,mpirun}`; `./uw build` after source changes. diff --git a/.claude/skills/adaptive-meshing/SKILL.md b/.claude/skills/adaptive-meshing/SKILL.md new file mode 100644 index 00000000..429efbb3 --- /dev/null +++ b/.claude/skills/adaptive-meshing/SKILL.md @@ -0,0 +1,361 @@ +--- +name: adaptive-meshing +description: The canonical, workable recipe for Underworld3 moving-mesh / adaptive-mesh convection (annulus stagnant-lid, faults, free surface). Reach for THIS first when setting up any model with a deforming or adapted mesh — it encodes the combination that does not blow up, tangle, or inject spurious energy, and explains the failure modes so you don't re-derive them. Use before choosing movers, free-slip BCs, restart, or field-transfer options. +--- + +# adaptive-meshing + +The one workable combination for UW3 moving/adaptive-mesh convection, distilled +from many sessions that each re-picked options and stomped on each other's +defaults. **Start from this recipe; change one thing at a time and verify.** + +Reference implementation (current, validated): the **`underworld3.workflows` +adaptive-convection example** — +`docs/examples/workflows/adaptive_convection/` on the `feature/adaptive-convection` +worktree (`config.py`+`simulate.py` no-fault; `fault_config.py`+`fault_simulate.py` +fault; `diagnostics.py`, `render.py`, `compare.py`). Express adaptive runs as a +WORKFLOW (a `WorkflowConfig` + `@workflow_step` DAG + `Run`), NOT a monolithic +driver. The older `scripts/fault_convection_adapt_loop.py` (feature/fault-convection) +is superseded — its ideas are folded into the workflow + this skill. +Companion: the `uw-visualisation` skill for rendering results. + +**Choosing the paradigm:** THIS skill is the **mover** (node movement / +equidistribution, `smooth_mesh_interior`) — the mesh deforms to follow a field. For +**local refinement** instead (`mesh.adapt(engine="nvb")` returns a refined CHILD; a +fault resolved by a fine band + custom-P FMG + rotated free-slip + dynamic topography ++ advection-diffusion), use the **`adapt-on-top-faults`** skill. Different tools — +don't mix them. + +--- + +## Mover quick-start (copy-paste — this is the hard-to-discover bit) + +The mesh mover is `uw.meshing.smooth_mesh_interior`. Minimal correct setup to +adapt a mesh to a field `T` each step: + +```python +import underworld3 as uw + +# metric from |grad T|: refinement = finest:coarsest cell-size ratio (~5). +# Use refinement=R, NOT strategy= (strategy caps at ~2 and under-grades). +rho = uw.meshing.metric_density_from_gradient( + mesh, T, refinement=5, coarsening="auto", metric_choice="front-following") + +# move the mesh — mmpde: variational, non-folding, clusters AND aligns cells. +# It OWNS field transfer (remaps T + SLCN history, fires on_remesh hooks). +uw.meshing.smooth_mesh_interior( + mesh, metric=rho, method="mmpde", + method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), # mmpde's OWN kwargs + slip_surfaces=True, # boundary nodes slide tangentially (parallel-safe) + skip_threshold=0.9) # skip the move when the mesh is already aligned +``` + +For a **sharp feature (fault)** pass an anisotropic SPD TENSOR metric instead of +the scalar `rho` (thin ACROSS the feature normal n) and bake a gmsh base: + +```python +import sympy +n = sympy.Matrix([nx, ny]) # constant fault-normal unit vector +d = dfac.sym[0] # DIRECT unsigned distance field (P1) +M = rho * sympy.eye(2) + (Rf**2 - 1.0) * sympy.exp(-(d/w)**2) * (n * n.T) +# mesh built with: uw.meshing.Annulus(..., refine_lines=[xy], refine_size_min=smin) +uw.meshing.smooth_mesh_interior(mesh, metric=M, method="mmpde", + method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), + slip_surfaces=True, skip_threshold=None) # tensor metric: do the skip check yourself +``` + +Pitfalls that make it "not work": `method="anisotropic"`/`"ot"` (shred/sliver — use +mmpde); injecting `relax`/`n_outer` (starves mmpde's CG); `strategy=` instead of +`refinement=R` (under-grades); a scalar bump for a fault (refines a fat corridor, +leaves the centre coarse); signed `Surface.distance.sym` for `d` (bleeds along the +line extension — use a direct unsigned distance). Full rationale + the rest of the +recipe (BCs, restart, field transfer, cadence) below. + +--- + +## The canonical recipe (defaults that work) + +### 1. Mover — mmpde +`uw.meshing.smooth_mesh_interior(mesh, metric=..., method="mmpde", +method_kwargs=dict(step_frac=0.2, accel="cg", momentum=0.0), slip_surfaces=True)`. +- mmpde = Huang–Kamenski variational, **non-folding** (energy → ∞ as detJ → 0), + clusters AND aligns cells. It is the only clean mover. +- **NOT** the `anisotropic` mover (shreds/freezes on static features) or `OT` + (slivers — OT is optimal transport, sliver-prone, not a route around the cap). +- Do **not** inject `relax`/`n_outer` into mmpde (those are the anisotropic + mover's knobs and starve mmpde's internal CG). + +### 2. Metric +- Thermal: `metric_density_from_gradient(mesh, T, refinement=R, + metric_choice="front-following")` — `refinement=R` (≈5) is the finest:coarsest + grading ratio; named `strategy=` caps at ~2 and under-grades. R≈5 extracts ~all + the grading the node budget/layout allows; don't over-tune R (benign no-op above + budget). +- Fault / sharp feature: a **hand-built anisotropic SPD tensor** + `M = ρ·I + (Rf²−1)·exp(−(d/w)²)·n nᵀ` (thin ACROSS the feature normal n). A + scalar bump refines a fat isotropic corridor and leaves the centre-line coarse. +- Use a DIRECT unsigned distance field (geometry_tools) for `d`, NOT the Surface's + signed `.distance.sym` (its zero-contour bleeds the metric along the line + extension). + +### 3. Creation vs maintenance (the cap) — and the gmsh base COMPOUNDS +mmpde **cannot create** strong refinement from a uniform mesh — it saturates at a +fixed-topology cap (~1.8× on the fault, re-measured), because the mover has a FIXED +node budget (it redistributes, never adds nodes). To go finer you need MORE NODES, +which only gmsh can add at construction. **Bake refinement into the gmsh base** +(`Annulus(refine_lines=[xy], refine_size_min=...)` — now real, see the Faults +section) and the mover doesn't just MAINTAIN it: the extra gmsh nodes lift it off +its budget cap so it **compounds** (gmsh f2 base 0.44 → mover 0.29; f3 base 0.30 → +mover 0.19 ≈ 5× finer). Measured (Ra1e6 Rf8 res24, all folded=0): +uniform 0.55 (~1.8×) → gmsh-f2 0.29 (~3.4×) → gmsh-f3 0.19 (~5×). Judge by +fault/bulk nearest-neighbour spacing RATIO, never by global misalignment. + +### 4. Cadence — adapt as often as you like (forced every step is FINE) +Adapt every step or every few — on a CORRECT build (see §5) the mover converges +dead-flat and forced every-step adaptation under vigorous convection is stable +(validated: 39/39 forced adapts, mesh folded=0, area-ratio ~14 flat). Skipping +when aligned just saves cost. (An earlier claim that forcing adaptation +"tangles / over-injects energy" was WRONG — that was the deformed-eval bug in §5.) + +### 5. THE bug that wrecked adaptation (holes) — and the fix +SYMPTOM: giant empty cells / holes in the adapted mesh, intermittent area-ratio +spikes, convection wrecked. ROOT CAUSE (proven): `uw.function.evaluate` +**mis-locates points on a deformed mesh** (the nav kd-tree `mesh._nav_coords` was +captured from the ORIGINAL coords and never refreshed), so the metric — built with +strictly-positive nodal values — evaluates to **NEGATIVE garbage** (even at its own +DOFs) → non-SPD → the mover wrecks the mesh. FIX (both): +- **`8a9d2ff2`** (refresh `_nav_coords` + projected normals on every deform) — + the real fix; makes `function.evaluate`/`points_in_domain` track deformation. +- **Monotone RBF metric bake** (in `_winslow_mmpde`): Shepard-interpolate the + metric from its **positive nodal values** — a convex average is guaranteed ≥0 + (monotone) + fast (no cell-location). Use RBF for the metric; it doesn't need + high-precision eval. +With these, the metric stays positive/SPD; #259's SPD-floor never fires (harmless). +The mover, `accel`, and `refinement=R` (e.g. R=5) are all FINE — they were +red-herring symptoms of the eval bug. Always confirm a CLEAN BUILD first +(`./uw build`; `md5 site-packages/.../smoothing.py == src`); a stale build is its +own cause of holes. + +### 6. Stokes free-slip +Penalty (`add_natural_bc(KFS·v·n·n)`, KFS≈1e6) or Nitsche γ=10 both work on a +UNIFORM mesh. Nitsche preferred for sharp fault corners. Do NOT diagnose free-slip +with nodal v·n: Nitsche enforces v·n=0 weakly, so nodal v·n is large even when +correct — use vrms (kinetic energy). (An earlier "warm-restart × Nitsche → blow-up, +use cold restart" diagnosis was WRONG / confounded by the §5 eval bug + a stale build.) + +**On an ADAPTIVE mesh, Nitsche free-slip γ=10 is TOO SOFT → intermittent vrms +spikes at adapt steps** (under-enforcement, NOT a blow-up: single-step velocity +garbage that T survives because the solve is cold-restarted — e.g. vrms +144→8887→144). The mover pins the boundary and refines just beneath it, creating +high-aspect-ratio near-boundary cells whose Nitsche inverse-estimate constant needs +a bigger penalty. Two fixes, BOTH good (corroborated across sessions): +- **Nitsche γ=100** (not 10) for a free-slip top on an adaptive mesh — clean, smooth + signal. This is INDEPENDENT of the penalty's mesh-size scaling: local vs global `h` + are equivalent here (both garbage at γ=10, both clean at γ=100). +- **Don't fully pin the boundary in the mover** — let it tangential-slip + (`slip_surfaces=True`, the §1 default), NOT `pinned_labels=[Upper,...]`. Avoiding + the distorted near-boundary cells lets γ=10 work. Pin a boundary only when you must + hold a prescribed shape (e.g. a free-surface height the integrator just set). + +Aside — free-SURFACE held-lid is the OPPOSITE problem: the GLOBAL `h = +get_min_radius` drifts BELOW the surface cell as the interior refines → the Nitsche +penalty OVER-stiffens → a spurious one-step surface "mountain" (vhmax spike). Fix = +a LOCAL per-cell `h` via `mesh.cell_size()` (deformation/adaptation-tracking); +`add_nitsche_bc(..., local_h=True)` is the default (PR #275). So held-lid wants LESS +penalty (local-h), free-slip-adaptive wants MORE (γ=100) — different knobs. + +### 7. Advection + timestep +`AdvDiffusionSLCN(mesh, u_Field=T, V_fn=V.sym, theta=1.0, monotone_mode='clamp')`. +theta=1.0 (backward Euler) for stability; monotone clamp bounds SL overshoot; +`V_fn=V.sym` is the PHYSICAL velocity. **dt can be LARGE — SLCN is unconditionally +stable; the smallest/median adapted cell does NOT have to govern dt.** Use +`estimate_dt(percentile=50)` × a multiplier (e.g. dt_mult 3–5) to advance physical +time faster (needed to develop convection in stiff stagnant-lid regimes). + +### 8. Rheology +- **isotropic** (linear) ⇒ `snes_type=ksponly` (one exact KSP solve; default + newtonls rejects steps on vigorous flow). Cheap. Use for resolved features. +- **TI / anisotropic** weak fault (`TransverseIsotropicFlowModel`: + `shear_viscosity_0=η_FK`, `shear_viscosity_1=η_weak`, `director`=fault normal): + keep the **default newtonls** (NOT ksponly — ksponly converges to the WRONG + answer); use **penalty** free-slip (Nitsche trips the anisotropic-Jacobian bug); + **GAMG** (FMG ~7× slower). Measured on the gmsh-resolved + one-sided base + (Ra1e6 Δη1e3): ran clean to t=0.06, folded=0, vrms→20, **~10×** the isotropic + cost per step (GAMG eats the 1000× anisotropic contrast — well under the feared + 20×). The TI fault visibly steers the flow (persistent recirculation at the trace). +- Viscosity-bearing fields are **P0/P1 only** (positivity; higher order overshoots). + FK viscosity `η = exp(θ(1−T))`, θ = ln(Δη). Floor a weak zone, don't multiply → 0. + +### 9. Build discipline +- `./uw build` after ANY source change; verify `uw.__file__` is in site-packages + and (when debugging) that `site-packages/.../smoothing.py` matches `src/...`. A + **stale mover build is a top cause of "giant empty elements / holes"** in the + adapted mesh. +- **NEVER** `pip install -e .` (contaminates all envs). Run from inside the worktree + (`pixi run -e amr-dev`). + +--- + +## Faults — the gmsh-resolved, on-fault, one-sided recipe (validated 2026-06-23) + +The full fault recipe, hard-won. Reference implementation: the +`underworld3.workflows` adaptive-convection example +(`docs/examples/workflows/adaptive_convection/{config,fault_config}.py`, on the +`feature/adaptive-convection` worktree — built on the FIXED mover base). Use the +WORKFLOW system, not a monolithic driver. + +### gmsh line refinement (`refine_lines`) — NOW IMPLEMENTED in `Annulus` +```python +uw.meshing.Annulus(radiusOuter=1, radiusInner=0.5, cellSize=1/24, + refine_lines=[xy], # list of (N,2) polylines (model coords) + refine_size_min=cellSize/3, # cell size ON the line (factor 2-3 is plenty) + refine_dist_min=0.02, refine_dist_max=0.12) # size ramps back to cellSize +``` +A gmsh **Distance + Threshold** field along the polyline; INTERIOR points are +embedded so nodes land ON the line. Backward-compatible (default `None`). It is a +core meshing change → lands as its own small meshing PR, separate from the workflow +example. (Before this session it was only *called* in old scripts and never existed +— don't trust `refine_lines` on any branch but the one carrying this commit.) + +### Keeping refinement ON the fault (the metric-composition trap) +The fault metric is `M = ρ·I + (Rf²−1)·exp(−(d/wₙ)²)·n nᵀ` with the isotropic SIZE +density `ρ`. **Do NOT fuse the fault density into ρ_T by PRODUCT** — the cold +surface thermal BL (ρ_T ~ R^d ~ 20–25) out-competes the fault near its top and the +refinement drifts ABOVE the fault, starving the deep fault ("seems to repel"): +- Use `ρ = max(ρ_T, fault_ρ)` (NOT `ρ_T · fault_ρ`). +- Make `fault_ρ = 1 + amp·gauss` with **amp > ρ_T** (~25 at R=5) so the fault wins + the max along its whole length. +- DIAGNOSE this by comparing the gmsh BASE (step 0, refinement centered on the + fault — correct) vs the DEVELOPED mesh (drifted above) → it's the MOVER's metric, + not gmsh. Render the mesh with the **fault trace overlaid** (`render.py --fault`). + +### One-sided fault influence (the clean control) +Even with max+amp, the symmetric metric DEMANDS both flanks while realized nodes +drift to the hanging wall. Make it one-sided: +- Store the **SIGNED** distance in `dfac` (the gaussians square it, so magnitude is + unchanged — the sign only feeds a `0.5(1+tanh(m·d/w))` gate). Probe the + radially-outward side once to define "upper" regardless of the distance tool's + orientation convention. +- `fault_metric_side` (both/upper/lower) gates the refinement; `fault_rheology_side` + gates the weak zone. **`both=upper`** is the physical recipe: a one-sided + hanging-wall damage zone with the mesh refined on the same side (refinement and + rheology coincide; gmsh-f3 → fault/bulk ~0.19, folded=0). `metric_side=lower` + instead pulls refinement onto the footwall to counter the upward drift. + +### The wedge fill (anti-collision) +The fault pull and the surface-BL pull compete for the coarse cells in the radial +sliver BETWEEN them. `fault_wedge=True` gmsh-fills that wedge (sample radial +segments from each fault point up to the surface, add as a second `refine_lines` +point set) so both pulls have their own budget and merge into one coherent fine +wedge instead of colliding. + +### Weak zone (rheology) — geometric blend + gaussian PEAKED on the fault (2026-06-24) +The weak-zone viscosity blends `η_FK` (background) and `floor` (fault) by the +influence `f`∈[0,1] (P1, positive). Get it right with TWO rules (verify by +reconstructing + rendering the REALIZED η field — `render_fields.py` — and the +combined `η_FK(T)^(1−f)` field; never assume the floor is reached): +- **GEOMETRIC blend `η_weak = η_FK^(1−f)·floor^f`** (NOT arithmetic + `η_FK·(1−f)+floor·f`, whose `(1−f)` term leaks the stiff-lid background through + → η≈8 at f=0.97 in a 1000× lid). Geometric reaches the floor genuinely (η≈1.2 at + f=0.97) AND ties the contrast to the LOCAL η_FK — so the fault automatically + bites hardest where it cuts the cold stiff lid (physically correct), nothing in + the hot interior. +- **`f` must be PEAKED on the fault (gaussian), NOT a TOPHAT block.** A top-hat + makes a uniform weak BLOCK; its sharp edges mean strain follows ∇f (TWO parallel + lines at the block edges, not on the fault), and a one-sided block is offset to + the hanging wall (η_1=1 core sits ABOVE the drawn line). Use a **gaussian + `f=exp(−(d/w)²)`** (peak f=1 ON the fault → η_1=1 on the drawn line, strain + localizes INTO the slot as a single feature). side=both = symmetric; side=upper = + hanging-wall halo (gaussian taper up, sharp footwall recovery — NO halving gate). + Centre the METRIC too (`metric_side=both`) so nodes refine on the line. +THERMAL CONTRAST (Louis's insight, confirmed): even a symmetric gaussian gives a +much bigger η-contrast on the COLD (upper/surface) side than the warm (footwall) +side, because η_FK rises ~1000× toward the surface — so the fault's dynamical +prominence naturally concentrates in the cold lid. This is a feature, not a bug. +Verified (gmsh-f3, gaussian width=0.025, side+metric=both): weak zone centred on +the drawn fault on the adapted mesh, folded=0. TUNE AT STEP 0 (build mesh+fields, +no solve — fast; plot the 1D η_1(d) profile + the 2D field). COST: genuine 1000× +TI contrast ~25–145 s/step (cold-start steps slow, ~25 s once developed; GAMG). +For TI see §8. (`fault_config.py`: `fault_profile=gaussian` (default still tophat — +pass gaussian), geometric blend in `create_solvers`; `render_fields.py` light maps.) + +--- + +## Failure modes — symptom → cause → fix + +| Symptom | Cause | Fix | +|---|---|---| +| Giant empty cells / holes in adapted mesh; intermittent area-ratio spikes | `function.evaluate` mis-locates on deformed mesh → metric → negative/non-SPD (the real bug). OR stale build. | `8a9d2ff2` + monotone RBF metric bake; `./uw build` + verify md5 site-packages==src | +| Decays when it should convect | over-diffusion OR under-resolution OR dt too small to develop | check a resolved arbiter (finer space+time); raise node budget; **larger dt** (dt_mult 3–5) | +| Fault won't refine under convection | mmpde creation cap; field gives no signal in cold lid | gmsh `refine_lines` base — the gmsh nodes let mmpde COMPOUND past the cap (~5×) | +| Refinement drifts ABOVE the fault / "repels", deep fault starved | fault density fused by PRODUCT with ρ_T → thermal BL out-competes the fault near the surface | `ρ = max(ρ_T, fault_ρ)`; `amp > ρ_T` (~25); or one-sided `fault_metric_side`; +`fault_wedge` | +| Refinement on both flanks but you want one side | symmetric (unsigned-distance) metric | signed `dfac` + `fault_metric_side`/`fault_rheology_side` tanh gate (`both=upper` = physical) | +| TI fault solve: ksponly gives wrong answer | ksponly skips the Picard the inexact GAMG inner solve needs | keep default newtonls; penalty free-slip; GAMG (~10× isotropic cost, fine) | +| free-slip ADAPTIVE: vrms spikes at adapt steps (e.g. 144→8887→144), T stays bounded | Nitsche γ=10 too soft on the distorted near-boundary cells the pinned-top+interior-refine creates (under-enforcement) | **Nitsche γ=100** (not 10); OR don't pin the boundary (tangential-slip `slip_surfaces=True`). h-scaling (local vs global) is moot here | +| free-SURFACE held-lid: spurious one-step surface "mountain" / vhmax spike | global `h=get_min_radius` drifts below the surface cell as interior refines → Nitsche OVER-stiff | LOCAL per-cell h: `add_nitsche_bc(local_h=True)` (default, PR #275) = `mesh.cell_size()` | + +**Diagnose by:** vrms (KE), Nu (`BdIntegral` surface flux), fault/bulk NN-spacing +RATIO (cKDTree), folded-element count + min cell area. Compare runs **at matched +physical time t** (dt differs between meshes), not by step number. When unsure +whether behaviour is physical, build a **resolved arbiter** (uniform mesh finer in +BOTH space and time) — if it agrees with one candidate, that's the truth. + +## Diagnostics (in the workflow example, reusable) +`diagnostics.py`: `mesh_quality` (folded / area-ratio / aspect), `nn_spacing_ratios` +(BL + fault/bulk), `NusseltSurface`, `vrms`, `History`. `render.py`: T+mesh+ +streamlines on the `Run` layout, **`--fault`** overlays the fault trace (read from +the run manifest) + **`--focus-fault`** auto-crops on it + **`--mesh-only`** for the +clean mesh, **`--all`** for every frame. `compare.py`/`fault_refine_plot.py`: +matched-physical-time comparison + fault/bulk-ratio time series. +**Rendering long runs as they go:** a completion-only Monitor is NOT enough — arm a +Monitor that POLLS for new `run.mesh.NNNNN.xdmf` checkpoints and emits the index, so +you render each step as it lands. +**Checkpoint an ADAPTIVE mesh with `meshUpdates=True`** (per-step geometry) or the +saved frames pair deformed fields with stale step-0 geometry. + +## The verified canonical command + +Requires the fix (§5): `8a9d2ff2` (deformed-mesh point-location) + the monotone RBF +metric bake in `smoothing.py`. Validated 2026-06-22 at vigorous Ra1e6/Δη1e3 +stagnant-lid: forced adapt EVERY step (39/39), larger dt, → vrms→18, Nu→1.86, +|v|max 61, mesh CLEAN every frame (folded=0, area-ratio ~14 flat), no abort. + +```bash +# no-fault baseline (the workflow CLI; one flag per config field) +pixi run -e amr-dev python docs/examples/workflows/adaptive_convection/simulate.py \ + --output-dir ~/+Simulations//baseline \ + --rayleigh 1e6 --delta-eta 1e3 --cellsize 0.0417 \ + --resolution-ratio 5 --adapt-every 1 --dt-mult 4 --max-steps 80 --max-t 0.06 + +# resolved fault: gmsh base (factor 3) + on-fault + one-sided hanging wall + TI +pixi run -e amr-dev python docs/examples/workflows/adaptive_convection/fault_simulate.py \ + --output-dir ~/+Simulations//fault \ + --rayleigh 1e6 --delta-eta 1e3 --cellsize 0.0417 --resolution-ratio 5 \ + --fault-base-smin 0.0139 --fault-anisotropy 8 \ + --metric-combine max --fault-refine-amp 25 \ + --fault-rheology-side upper --fault-metric-side upper \ + --rheology ti --dt-mult 4 --max-steps 40 --max-t 0.06 +``` + +Key choices, verified: +- **`--resolution-ratio 5`** (R=5) is fine — the mover handles it on a correct build. +- **`--adapt-every 1`** force adapt every step (the strongest mesh test). +- **`--dt-mult 4`** — larger dt for STABILITY is fine (SLCN unconditional); but it + costs transient ACCURACY (over-diffusive backward-Euler DELAYS the convective + onset vs a resolved arbiter — dt×1.5 recovers it). Use small dt_mult for faithful + transients, large to reach quasi-steady fast. +- **`--freeslip penalty`** (default) — REQUIRED for TI (Nitsche trips the + anisotropic-Jacobian bug). It's the raw velocity penalty `kfs·(v·n)n`. +- Fault: `--fault-base-smin` (gmsh resolve), `--metric-combine max` + + `--fault-refine-amp 25` (keep refinement ON the fault), `--fault-*-side` (one-sided), + `--rheology ti` (real fault). Drop the fault flags for the no-fault control. + +Render with `render.py` (`--fault --focus-fault` to see the trace + refinement +coincide; `--mesh-only`; `--all`). Judge the mesh by folded/area-ratio, the physics +by vrms/Nu, ALWAYS at matched physical time. + +## Related memory +`project_adaptive_convection_as_workflow` (THIS session: workflow port, gmsh +refine_lines, on-fault/one-sided/wedge, TI, dt-accuracy), `project_mmpde_holes_real_root_cause`, +`project_fault_refine_fixed_topology_cap`, `project_uw_workflow_landing`, +`feedback_debug_adaptive_solver_method`, `project_fault_convection_working_settings`. diff --git a/.claude/skills/cetz-figures/SKILL.md b/.claude/skills/cetz-figures/SKILL.md new file mode 100644 index 00000000..976f4c97 --- /dev/null +++ b/.claude/skills/cetz-figures/SKILL.md @@ -0,0 +1,138 @@ +--- +name: cetz-figures +description: Build schematic / labelled-geometry figures for underworld3 papers using Typst + cetz. Use when the figure is primarily about topology, annotation, and math-typeset labels (meshes, solver diagrams, flow charts). Prefer Python → SVG → `#image()` for data-heavy figures (fields, colormaps, arrow plots) instead. +--- + +# cetz-figures + +Scaffold for Typst/cetz figures in `publications/**/figures/` alongside the +existing `arrays-sync-flow.typ`. This skill exists because upstream Claude +sessions draft cetz blind — this one actually compiles. + +## When to use cetz + +- Mesh schematics with a few labelled triangles / vertices / control points. +- Solver / data-flow diagrams (see `arrays-sync-flow.typ` in this repo). +- Anything where labels should render in the paper's math/text fonts. +- Anything that benefits from recompiling with the paper. + +## When to use something else + +- **Data-heavy plots** (scalar fields, colormaps, quiver plots, anything with + dense per-pixel or per-cell data) — generate SVG from Python/matplotlib, + include via `#image("foo.svg")`. cetz will fight you here. +- **Geometry computation** (Delaunay, intersections, interpolation) — do it in + Python offline, emit JSON with the shape `{"vertices": [...], + "triangles": [...], ...}`, let Typst just draw. + +## When NOT to use TikZ + +Evaluated and rejected: +- Slower compile than Typst. +- Drags in a LaTeX toolchain that isn't otherwise required by the project. +- No Typst math-font advantage over cetz for labels in our paper context. + +Keep TikZ in back pocket only if a co-author insists on TikZ source. + +## Before you draw (thesis-first discipline) + +When a user hands you a figure request — especially one replacing an +existing ASCII sketch, whiteboard photo, or reference figure — **do not +start transcribing**. The source artefact is a hypothesis about what +to communicate, not a specification. Before opening cetz: + +1. **State the thesis in one sentence.** What is the figure arguing? + If you can't write it plainly, you don't yet understand the figure. + Ask the user to articulate it. + +2. **Honestly audit the source.** If the original ASCII / sketch is + being replaced, it's often replaced *because it doesn't land well*. + Say what's broken about it before proposing the replacement — the + user often agrees and the new figure can do more than the old. + +3. **Enumerate design decisions as explicit questions, not assumptions.** + For a curved-boundary normals figure that's typically: + - Geometry (circle / ellipse / arc span / zoom level) + - Sampling density (number of facets, quadrature points per facet) + - Overlay vs. side-by-side + - Whether to show error quantitatively (arcs, annotations) or leave + it as a visible angle + - Where the figure lives in the repo (which doc / which branch) + Present a proposed interpretation with the decisions flagged; let + the user resolve them before you compile. + +4. **Only then open cetz.** Iterate visually — the discipline above is + about not committing to a design prematurely, not about planning + exhaustively. Once you start, compile often. + +The user's phrase "I'm not quite sure what this is intended to +illustrate" is the canonical trigger for this discipline. If you hear +it (or catch yourself about to transcribe without checking), stop and +do the four steps above. + +## Key gotchas (hit during iteration — not hypothetical) + +1. **Don't name helper parameters `anchor`.** cetz internals treat `anchor:` + specially and you get a cryptic panic: `"Unknown anchor 'anchor' for + element 'none'"`. Use `align-to` or similar. See `cetz-cheatsheet.md`. + +2. **Clipping is a Typst concern, not a cetz one.** cetz has no `\clip`. + Wrap the canvas in `#box(clip: true, width: ..., height: ..., ...)` and + draw slightly oversized inside the canvas — the box clips the overflow. + +3. **Painter's algorithm — order matters.** Draw the background first, the + highlight second. No z-index exists. Verified in `mesh-demo.typ`. + +4. **Semi-transparency via `rgb(r, g, b, a)`** (alpha 0–255) or + `color.transparentize(col, 50%)`. Both work for fill and stroke. + +5. **Math in labels just works.** `content(pos, $v_1$)` renders in the + document math font. No escape hatch needed. This is a real cetz win over + SVG. + +6. **`import cetz.draw: *` inside the canvas closure.** Without it, `line`, + `circle`, `content` aren't in scope. Helper functions that draw need + their own `import cetz.draw: *` line inside. + +## Project layout pattern + +Each blog post or paper section gets its own subdirectory under +`figures/`, so a post's figures travel together: + +``` +publications/blog-posts/figures/ +└── / + ├── .typ # cetz drawing + ├── .png # committed output + ├── -data.json # (optional) precomputed geometry + └── generate-.py # (optional) Python that writes the JSON +``` + +Concrete example: `publications/blog-posts/figures/finding-particles/` +holds `mesh-demo.*` and `domain-demo.*` for the post +`finding-particles.md`. + +The JSON intermediate is the forward bridge to underworld3 — see +`underworld-bridge.md`. + +## Reference files + +- `cetz-cheatsheet.md` — what worked from memory vs. needed lookup. +- `underworld-bridge.md` — JSON schema for future `uw.meshing` export. +- `examples/` — self-contained copies (each with `.typ`, `.png`, `.json`, + and generator `.py`) of the figures this skill is scaffolded from. + These are snapshots; the live versions may have drifted if a post or + doc was iterated on further. + - `mesh-demo.*` — element-level point-in-cell test. + Live: `publications/blog-posts/figures/finding-particles/`. + - `domain-demo.*` — parallel domain centroid ambiguity. + Live: `publications/blog-posts/figures/finding-particles/`. + - `facet-vs-true-normals.*` — facet normal vs. smooth-surface normal + on a curved boundary. Live: + `docs/advanced/figures/curved-bc/`. + +## Canonical reference in the repo + +`publications/blog-posts/figures/arrays-sync-flow.typ` — prior cetz figure +in the repo, established version (0.3.4) and house style (hex colours, +helper-function pattern). Follow its conventions. diff --git a/.claude/skills/cetz-figures/cetz-cheatsheet.md b/.claude/skills/cetz-figures/cetz-cheatsheet.md new file mode 100644 index 00000000..4d3c190c --- /dev/null +++ b/.claude/skills/cetz-figures/cetz-cheatsheet.md @@ -0,0 +1,312 @@ +# cetz cheatsheet — honest edition + +Separates **what I used correctly from memory** (so you probably can too) +from **what actually needed diagnosis** (so you should check the manual +at https://cetz-package.github.io/docs/ rather than trust the pattern). + +Version: cetz 0.3.4, Typst 0.13.1. Don't assume these patterns hold for +later cetz versions — the API has been unstable between minor versions. + +**External reference worth checking first:** + — curated gallery of ~130 scientific +diagrams with source (~120 in cetz, ~110 in TikZ). Heaviest coverage is +theoretical physics and ML; lighter on FE/meshes specifically, but +strong patterns for coordinate frames, Euler angles, 3D orientation, +and lattice structures. Worth browsing before inventing a new figure +from scratch — the answer to "how do people draw X in cetz" is often +sitting there. + +## Worked from memory (verified correct) + +```typst +#import "@preview/cetz:0.3.4" + +#cetz.canvas(length: 4cm, { + import cetz.draw: * + + // Straight line / polyline / closed polygon + line((0, 0), (1, 0), (1, 1), close: true, + fill: rgb(70, 120, 180, 60), + stroke: 1.2pt + rgb("#1f3a6b")) + + // Filled dot + circle((0.5, 0.5), radius: 0.02, fill: black, stroke: none) + + // Label at a position — math content renders in document math font + content((0.5, 0.5), $v_1$, anchor: "west") +}) +``` + +- `#box(clip: true, width: 8cm, height: 8cm, cetz.canvas(...))` — the + *Typst* box does the clipping. cetz itself doesn't clip. +- `rgb("#b8b8b8")` for hex colours, `rgb(r, g, b, a)` for RGBA (alpha 0–255). +- `json("file.json")` from outside the canvas, then `.at(i)` and field + access inside. Values flow through into cetz coordinates cleanly. +- `for tri in mesh.triangles { line(...) }` — plain Typst for-loops work + inside the canvas closure. + +## Needed diagnosis (check the manual if you do this) + +### Parameter-name collision with `anchor` + +```typst +// BROKEN — panic: "Unknown anchor 'anchor' for element 'none'" +#let dot(p, label, anchor: "west") = { + circle(p, radius: 0.02, fill: black) + content(p, label, anchor: anchor) +} +``` + +```typst +// FIX — rename the parameter +#let dot(p, label, align-to: "west") = { + circle(p, radius: 0.02, fill: black) + content(p, label, anchor: align-to) +} +``` + +Diagnosis cost: one compile failure, ~5 minutes of confusion, then the +rename. The error message does not point at your parameter. + +### Anchor values for `content()` + +Confirmed to work: `"north"`, `"south"`, `"east"`, `"west"`, +`"north-east"`, `"south-east"`, `"north-west"`, `"south-west"`, `"center"`. + +Not tried in this session: named anchors of other elements (e.g. +`"my-rect.east"`), element-specific anchors beyond compass. Verify +against the cetz anchors docs if you need those. + +### Content wider than the Typst box needs `align()` to re-centre + +`#box(clip: true, width: W, height: H, cetz.canvas(...))` does not centre +the cetz canvas if the cetz bounding box is larger than `W × H` (e.g. a +mesh drawn oversize for clipping). Typst default-aligns the canvas to +the top-left, which shifts the interesting region off-frame. Wrap: + +```typst +#box(clip: true, width: W, height: H, + align(center + horizon, cetz.canvas(length: L, { ... }))) +``` + +Encountered when the background mesh extent grew past `[-1, 1]` — +first compile looked fine because content fit; the next compile didn't. + +### Dashed / dotted strokes in cetz + +Pass a dict for the stroke: + +```typst +stroke: (paint: rgb("#1f3a6b"), thickness: 0.5pt, dash: "dotted") +stroke: (paint: rgb("#059669"), thickness: 0.7pt, dash: "dashed") +``` + +Valid `dash` values include `"dotted"`, `"dashed"`, `"densely-dotted"`, +`"loosely-dotted"`, etc. (Typst stroke spec, not a cetz extension.) + +### Line-to-viewport clipping is your job + +cetz does not clip lines at an imaginary viewport. If you want to +extend the line through `(p, q)` out to the edge of a `[-V, V]²` box, +compute the `t` values at each axis crossing, take the inner `t-enter` +and `t-exit`, and evaluate the line there. See `extend-to-viewport` in +`examples/mesh-demo.typ`. + +### Draw helpers that take `cetz.draw` primitives + +Helper functions defined **outside** `cetz.canvas({...})` need their +own `import cetz.draw: *` inside the function body: + +```typst +#let dot(p, label, ...) = { + import cetz.draw: * // REQUIRED — cetz.draw symbols not in scope otherwise + circle(p, ...) + content(...) +} +``` + +This feels like it shouldn't be necessary but is (cetz 0.3.4). + +## Probably wrong / would need checking + +These are *guesses* from adjacent Typst/cetz knowledge, not verified here: + +- `plot` module usage (axes, line plots, scatter plots). Upstream sessions + were shaky on this; assume the API has drifted and check the cetz plot + docs before using. +- `decorations` module (brace, waveform, etc.). Not touched in this session. +- Named-anchor references across draw calls via `name:` — syntax is + `element-name.anchor-name` but context behaviour needs testing. +- Gradient fills — not used, not verified. +- Transformations (`translate`, `rotate`, `scale`) — likely work but not + exercised here. + +## Compile workflow + +```bash +# First run downloads cetz and its deps (small, ~90 KB total) +typst compile mesh-demo.typ mesh-demo.png + +# PDF is default if no output specified +typst compile mesh-demo.typ # → mesh-demo.pdf + +# Watch mode for iteration +typst watch mesh-demo.typ mesh-demo.png +``` + +## Data-driven figures with JSON + +The preferred pattern for mesh/particle figures is to generate data in +Python (scipy.spatial.Delaunay, numpy), export as JSON, and render in +Typst/cetz. This separates the geometry computation from the drawing. + +```typst +#let mesh = json("mesh-data.json") +#let pt(i) = { + let v = mesh.vertices.at(i) + (v.at(0), v.at(1)) +} +// Loop over triangles +for tri in mesh.triangles { + line(pt(tri.at(0)), pt(tri.at(1)), pt(tri.at(2)), + close: true, fill: fill-col, stroke: stroke-mesh) +} +``` + +### Python generator pattern + +```python +from scipy.spatial import Delaunay +# Jittered equilateral lattice → Delaunay → JSON +tri = Delaunay(points) +data = { + "vertices": [[float(x), float(y)] for x, y in points], + "triangles": [[int(i), int(j), int(k)] for i, j, k in tri.simplices], + # ... domain membership, centroids, test points, etc. +} +json.dump(data, open("mesh-data.json", "w"), indent=2) +``` + +### Enumeration with index + +```typst +// .enumerate() gives (index, value) pairs — used for domain colouring +for (idx, tri) in mesh.triangles.enumerate() { + let d = get-domain(idx) + let fill-col = if d == "A" { fill-a } else { fill-b } + line(pt(tri.at(0)), pt(tri.at(1)), pt(tri.at(2)), + close: true, fill: fill-col, stroke: stroke-mesh) +} +``` + +### High-DPI PNG output + +```bash +typst compile figure.typ figure.png --ppi 300 +``` + +## Lessons from domain-demo figure + +- **Multiple domain colours**: Use low-alpha fills (25–45) so mesh edges + show through. Peripheral domains use even lower alpha to stay in background. +- **No heavy boundary lines needed**: Domain colours communicate the + partition without explicit boundary drawing. This looks more natural. +- **Smaller dots and tight labels** for a cleaner look: `radius: 0.028–0.032`, + label offsets of `(0.06, 0.06)` or `(0.08, 0.06)`. +- **Finer mesh** (SPACING ≈ 0.28, ~400 vertices, ~800 triangles) reads as + a realistic FE mesh. Coarser meshes (SPACING ≈ 1.25) work for + element-level diagrams. +- **View clipping**: Use `VIEW` parameter in JSON and `in-view()` check + to clip triangles at the edge of the rendered region. + +## Visual idioms that worked + +### Dashed-on-top-of-solid reveals coincidence + +When two arrows (or two paths) represent quantities that *sometimes* +coincide and *sometimes* diverge, draw the **dashed** one LAST so it +layers on top of the **solid** one. Where they overlap, the dashes +interleave with the solid underneath and both colours remain visible +— reading unmistakably as "they coincide here". If the solid were on +top, it would simply cover the dashed and the coincidence would be +invisible. + +Applied in `facet-vs-true-normals.typ`: the rust `mesh.Gamma` arrow is +drawn first, the green dashed true normal on top. At each facet +midpoint the green dashes sit over the rust — the coincidence is +obvious. At off-centre quadrature points both arrows are visible as +distinct directions. + +### Quadrature choice encodes pedagogy + +If your figure illustrates *within-element variation* of a quantity, +the quadrature rule you pick determines whether the variation is +visible at all: + +- **2-point Gauss–Legendre** (±1/√3 on [-1, 1]): both samples off + the element midpoint. Can show *average* per-element behaviour + but hides any quantity that's zero at the midpoint. +- **3-point Gauss–Legendre** (0, ±√(3/5)): one sample AT the + midpoint plus two off-centre. Reveals both the midpoint case and + the at-the-edge case in the same picture. + +Switching from 2-point to 3-point was the key move that made the +curved-boundary figure communicate "the error is zero at facet +midpoints" — not just "there is an error". + +### Legend in empty interior space beats per-arrow labels + +Per-arrow labels near the canvas edge fight cetz's bbox calculation +and often clip at the PNG boundary. If the figure has an obvious +empty interior region, put a small 2-row legend there: + +```typst +// Row 1 +line((lx, ly), (lx + 0.28, ly), stroke: gamma-stroke, + mark: (end: ">", fill: gamma-colour)) +content((lx + 0.36, ly), text(fill: gamma-colour, $hat(n)_Gamma$ + [ (label)]), + anchor: "west") +// Row 2, below +``` + +The short sample stroke uses the *same* stroke dict as the arrows in +the figure, so it looks identical. No leader lines, no per-arrow +crowding. + +### Short arrows read better than long + +Arrow lengths around **0.25–0.3 world units** (with canvas `length: +3cm`) stay visually subordinate to the geometry they annotate. +Longer arrows (0.4+) tend to dominate and obscure the mesh. Start +short; lengthen only if the arrow directions are genuinely hard to +read. + +### Explicit page dimensions + aligned canvas = reliable padding + +For docs figures that must look crisp in a PNG of specific dimensions: + +```typst +#set page(width: 12cm, height: 8cm, margin: 6pt) +#align(center + horizon, cetz.canvas(length: 3cm, { ... })) +``` + +This decouples the *rendered size* from cetz's (sometimes generous, +sometimes flaky) auto-bbox. Your content scales predictably and text +labels near the edges have room. Use it any time the figure goes +into a doc as a fixed-size image. + +## House style in this repo + +Figures live in per-post subdirectories under `publications/blog-posts/figures/`: +- `finding-particles/mesh-demo.typ` — element-level inside/outside test with control points +- `finding-particles/domain-demo.typ` — multi-domain parallel mesh with centroid ambiguity +- `arrays-sync-flow.typ` — flow diagram with nodes, regions, bezier connectors + (older figure, not yet moved into a per-post subdirectory) + +Patterns: +- Custom helper functions declared at the top of the canvas closure + (`node`, `region`, `elabel`, `dot`). +- Hex colours via `rgb("#...")`. +- Arrows via `mark: (end: ">", fill: black)`. +- Bezier connectors via `bezier(start, end, ctrl1, ctrl2, ...)`. +- Data from JSON files generated by companion Python scripts. diff --git a/.claude/skills/cetz-figures/examples/curved-bc-data.json b/.claude/skills/cetz-figures/examples/curved-bc-data.json new file mode 100644 index 00000000..0abb0757 --- /dev/null +++ b/.claude/skills/cetz-figures/examples/curved-bc-data.json @@ -0,0 +1,652 @@ +{ + "centre": [ + 0.0, + 0.0 + ], + "radius": 1.5, + "arc_points": [ + [ + 1.299, + 0.75 + ], + [ + 1.2858, + 0.7726 + ], + [ + 1.2721, + 0.7949 + ], + [ + 1.258, + 0.817 + ], + [ + 1.2436, + 0.8388 + ], + [ + 1.2287, + 0.8604 + ], + [ + 1.2135, + 0.8817 + ], + [ + 1.198, + 0.9027 + ], + [ + 1.182, + 0.9235 + ], + [ + 1.1657, + 0.944 + ], + [ + 1.1491, + 0.9642 + ], + [ + 1.1321, + 0.9841 + ], + [ + 1.1147, + 1.0037 + ], + [ + 1.097, + 1.023 + ], + [ + 1.079, + 1.042 + ], + [ + 1.0607, + 1.0607 + ], + [ + 1.042, + 1.079 + ], + [ + 1.023, + 1.097 + ], + [ + 1.0037, + 1.1147 + ], + [ + 0.9841, + 1.1321 + ], + [ + 0.9642, + 1.1491 + ], + [ + 0.944, + 1.1657 + ], + [ + 0.9235, + 1.182 + ], + [ + 0.9027, + 1.198 + ], + [ + 0.8817, + 1.2135 + ], + [ + 0.8604, + 1.2287 + ], + [ + 0.8388, + 1.2436 + ], + [ + 0.817, + 1.258 + ], + [ + 0.7949, + 1.2721 + ], + [ + 0.7726, + 1.2858 + ], + [ + 0.75, + 1.299 + ], + [ + 0.7272, + 1.3119 + ], + [ + 0.7042, + 1.3244 + ], + [ + 0.681, + 1.3365 + ], + [ + 0.6576, + 1.3482 + ], + [ + 0.6339, + 1.3595 + ], + [ + 0.6101, + 1.3703 + ], + [ + 0.5861, + 1.3808 + ], + [ + 0.5619, + 1.3908 + ], + [ + 0.5376, + 1.4004 + ], + [ + 0.513, + 1.4095 + ], + [ + 0.4884, + 1.4183 + ], + [ + 0.4635, + 1.4266 + ], + [ + 0.4386, + 1.4345 + ], + [ + 0.4135, + 1.4419 + ], + [ + 0.3882, + 1.4489 + ], + [ + 0.3629, + 1.4554 + ], + [ + 0.3374, + 1.4616 + ], + [ + 0.3119, + 1.4672 + ], + [ + 0.2862, + 1.4724 + ], + [ + 0.2605, + 1.4772 + ], + [ + 0.2347, + 1.4815 + ], + [ + 0.2088, + 1.4854 + ], + [ + 0.1828, + 1.4888 + ], + [ + 0.1568, + 1.4918 + ], + [ + 0.1307, + 1.4943 + ], + [ + 0.1046, + 1.4963 + ], + [ + 0.0785, + 1.4979 + ], + [ + 0.0523, + 1.4991 + ], + [ + 0.0262, + 1.4998 + ], + [ + 0.0, + 1.5 + ], + [ + -0.0262, + 1.4998 + ], + [ + -0.0523, + 1.4991 + ], + [ + -0.0785, + 1.4979 + ], + [ + -0.1046, + 1.4963 + ], + [ + -0.1307, + 1.4943 + ], + [ + -0.1568, + 1.4918 + ], + [ + -0.1828, + 1.4888 + ], + [ + -0.2088, + 1.4854 + ], + [ + -0.2347, + 1.4815 + ], + [ + -0.2605, + 1.4772 + ], + [ + -0.2862, + 1.4724 + ], + [ + -0.3119, + 1.4672 + ], + [ + -0.3374, + 1.4616 + ], + [ + -0.3629, + 1.4554 + ], + [ + -0.3882, + 1.4489 + ], + [ + -0.4135, + 1.4419 + ], + [ + -0.4386, + 1.4345 + ], + [ + -0.4635, + 1.4266 + ], + [ + -0.4884, + 1.4183 + ], + [ + -0.513, + 1.4095 + ], + [ + -0.5376, + 1.4004 + ], + [ + -0.5619, + 1.3908 + ], + [ + -0.5861, + 1.3808 + ], + [ + -0.6101, + 1.3703 + ], + [ + -0.6339, + 1.3595 + ], + [ + -0.6576, + 1.3482 + ], + [ + -0.681, + 1.3365 + ], + [ + -0.7042, + 1.3244 + ], + [ + -0.7272, + 1.3119 + ], + [ + -0.75, + 1.299 + ], + [ + -0.7726, + 1.2858 + ], + [ + -0.7949, + 1.2721 + ], + [ + -0.817, + 1.258 + ], + [ + -0.8388, + 1.2436 + ], + [ + -0.8604, + 1.2287 + ], + [ + -0.8817, + 1.2135 + ], + [ + -0.9027, + 1.198 + ], + [ + -0.9235, + 1.182 + ], + [ + -0.944, + 1.1657 + ], + [ + -0.9642, + 1.1491 + ], + [ + -0.9841, + 1.1321 + ], + [ + -1.0037, + 1.1147 + ], + [ + -1.023, + 1.097 + ], + [ + -1.042, + 1.079 + ], + [ + -1.0607, + 1.0607 + ], + [ + -1.079, + 1.042 + ], + [ + -1.097, + 1.023 + ], + [ + -1.1147, + 1.0037 + ], + [ + -1.1321, + 0.9841 + ], + [ + -1.1491, + 0.9642 + ], + [ + -1.1657, + 0.944 + ], + [ + -1.182, + 0.9235 + ], + [ + -1.198, + 0.9027 + ], + [ + -1.2135, + 0.8817 + ], + [ + -1.2287, + 0.8604 + ], + [ + -1.2436, + 0.8388 + ], + [ + -1.258, + 0.817 + ], + [ + -1.2721, + 0.7949 + ], + [ + -1.2858, + 0.7726 + ], + [ + -1.299, + 0.75 + ] + ], + "facet_vertices": [ + [ + 1.299, + 0.75 + ], + [ + 0.513, + 1.4095 + ], + [ + -0.513, + 1.4095 + ], + [ + -1.299, + 0.75 + ] + ], + "quadrature": [ + { + "pos": [ + 1.2105, + 0.8243 + ], + "facet_normal": [ + 0.6428, + 0.766 + ], + "true_normal": [ + 0.8265, + 0.5629 + ], + "facet_idx": 0 + }, + { + "pos": [ + 0.906, + 1.0798 + ], + "facet_normal": [ + 0.6428, + 0.766 + ], + "true_normal": [ + 0.6428, + 0.766 + ], + "facet_idx": 0 + }, + { + "pos": [ + 0.6016, + 1.3352 + ], + "facet_normal": [ + 0.6428, + 0.766 + ], + "true_normal": [ + 0.4108, + 0.9117 + ], + "facet_idx": 0 + }, + { + "pos": [ + 0.3974, + 1.4095 + ], + "facet_normal": [ + 0.0, + 1.0 + ], + "true_normal": [ + 0.2714, + 0.9625 + ], + "facet_idx": 1 + }, + { + "pos": [ + 0.0, + 1.4095 + ], + "facet_normal": [ + 0.0, + 1.0 + ], + "true_normal": [ + 0.0, + 1.0 + ], + "facet_idx": 1 + }, + { + "pos": [ + -0.3974, + 1.4095 + ], + "facet_normal": [ + 0.0, + 1.0 + ], + "true_normal": [ + -0.2714, + 0.9625 + ], + "facet_idx": 1 + }, + { + "pos": [ + -0.6016, + 1.3352 + ], + "facet_normal": [ + -0.6428, + 0.766 + ], + "true_normal": [ + -0.4108, + 0.9117 + ], + "facet_idx": 2 + }, + { + "pos": [ + -0.906, + 1.0798 + ], + "facet_normal": [ + -0.6428, + 0.766 + ], + "true_normal": [ + -0.6428, + 0.766 + ], + "facet_idx": 2 + }, + { + "pos": [ + -1.2105, + 0.8243 + ], + "facet_normal": [ + -0.6428, + 0.766 + ], + "true_normal": [ + -0.8265, + 0.5629 + ], + "facet_idx": 2 + } + ], + "radius_end": [ + 0.0, + 1.5 + ] +} \ No newline at end of file diff --git a/.claude/skills/cetz-figures/examples/domain-demo-data.json b/.claude/skills/cetz-figures/examples/domain-demo-data.json new file mode 100644 index 00000000..537920bb --- /dev/null +++ b/.claude/skills/cetz-figures/examples/domain-demo-data.json @@ -0,0 +1,6275 @@ +{ + "vertices": [ + [ + -2.4847, + -2.5034 + ], + [ + -2.1999, + -2.4889 + ], + [ + -1.9627, + -2.4734 + ], + [ + -1.6454, + -2.484 + ], + [ + -1.4008, + -2.5028 + ], + [ + -1.1072, + -2.4761 + ], + [ + -0.8119, + -2.4819 + ], + [ + -0.5432, + -2.5153 + ], + [ + -0.2569, + -2.5244 + ], + [ + 0.0383, + -2.4926 + ], + [ + 0.3145, + -2.5081 + ], + [ + 0.6064, + -2.478 + ], + [ + 0.8756, + -2.5171 + ], + [ + 1.1381, + -2.5255 + ], + [ + 1.4006, + -2.4897 + ], + [ + 1.7137, + -2.4738 + ], + [ + 1.9702, + -2.5073 + ], + [ + 2.2583, + -2.5174 + ], + [ + -2.3807, + -2.2589 + ], + [ + -2.0953, + -2.248 + ], + [ + -1.8035, + -2.2389 + ], + [ + -1.5088, + -2.268 + ], + [ + -1.2214, + -2.2404 + ], + [ + -0.9663, + -2.2694 + ], + [ + -0.6698, + -2.2777 + ], + [ + -0.4168, + -2.2851 + ], + [ + -0.1039, + -2.2483 + ], + [ + 0.1715, + -2.2418 + ], + [ + 0.4377, + -2.2537 + ], + [ + 0.6998, + -2.2791 + ], + [ + 1.0094, + -2.2591 + ], + [ + 1.2837, + -2.2427 + ], + [ + 1.5675, + -2.2545 + ], + [ + 1.8433, + -2.2685 + ], + [ + 2.0937, + -2.2611 + ], + [ + 2.384, + -2.2626 + ], + [ + -2.4802, + -2.0299 + ], + [ + -2.2447, + -2.0273 + ], + [ + -1.9516, + -2.006 + ], + [ + -1.6568, + -1.9991 + ], + [ + -1.3708, + -2.0203 + ], + [ + -1.0824, + -2.0337 + ], + [ + -0.8467, + -2.038 + ], + [ + -0.5275, + -2.0172 + ], + [ + -0.279, + -2.015 + ], + [ + 0.0005, + -2.004 + ], + [ + 0.297, + -2.0217 + ], + [ + 0.5689, + -2.0077 + ], + [ + 0.8523, + -2.0381 + ], + [ + 1.1186, + -1.9892 + ], + [ + 1.4429, + -2.0038 + ], + [ + 1.6869, + -1.9888 + ], + [ + 1.9956, + -2.0029 + ], + [ + 2.2572, + -2.0278 + ], + [ + -2.3826, + -1.75 + ], + [ + -2.0825, + -1.7892 + ], + [ + -1.8109, + -1.7681 + ], + [ + -1.5381, + -1.7526 + ], + [ + -1.2255, + -1.7602 + ], + [ + -0.9638, + -1.7654 + ], + [ + -0.6753, + -1.7641 + ], + [ + -0.4233, + -1.7773 + ], + [ + -0.1457, + -1.7729 + ], + [ + 0.1505, + -1.7924 + ], + [ + 0.4178, + -1.7676 + ], + [ + 0.7016, + -1.7487 + ], + [ + 1.0045, + -1.7811 + ], + [ + 1.2851, + -1.7993 + ], + [ + 1.5857, + -1.7735 + ], + [ + 1.8558, + -1.7959 + ], + [ + 2.1193, + -1.7731 + ], + [ + 2.4245, + -1.7685 + ], + [ + -2.5015, + -1.5431 + ], + [ + -2.2294, + -1.5289 + ], + [ + -1.9434, + -1.5568 + ], + [ + -1.6417, + -1.5079 + ], + [ + -1.4001, + -1.527 + ], + [ + -1.1219, + -1.5204 + ], + [ + -0.8323, + -1.5211 + ], + [ + -0.5273, + -1.515 + ], + [ + -0.282, + -1.5068 + ], + [ + 0.0049, + -1.556 + ], + [ + 0.3031, + -1.5373 + ], + [ + 0.5985, + -1.5128 + ], + [ + 0.8498, + -1.5047 + ], + [ + 1.1283, + -1.5292 + ], + [ + 1.4063, + -1.5056 + ], + [ + 1.6812, + -1.5555 + ], + [ + 1.9764, + -1.5025 + ], + [ + 2.2819, + -1.5161 + ], + [ + -2.3381, + -1.2655 + ], + [ + -2.0789, + -1.2979 + ], + [ + -1.7848, + -1.2785 + ], + [ + -1.5271, + -1.3103 + ], + [ + -1.2262, + -1.3009 + ], + [ + -0.9355, + -1.3021 + ], + [ + -0.7011, + -1.269 + ], + [ + -0.4194, + -1.3055 + ], + [ + -0.1144, + -1.2666 + ], + [ + 0.143, + -1.2982 + ], + [ + 0.4555, + -1.2611 + ], + [ + 0.72, + -1.3075 + ], + [ + 0.9728, + -1.3027 + ], + [ + 1.2594, + -1.2776 + ], + [ + 1.5388, + -1.2872 + ], + [ + 1.8509, + -1.283 + ], + [ + 2.1032, + -1.2705 + ], + [ + 2.4121, + -1.2742 + ], + [ + -2.5207, + -1.0661 + ], + [ + -2.1961, + -1.0508 + ], + [ + -1.9511, + -1.0457 + ], + [ + -1.6509, + -1.0196 + ], + [ + -1.392, + -1.0213 + ], + [ + -1.1266, + -1.042 + ], + [ + -0.8125, + -1.0671 + ], + [ + -0.5601, + -1.0496 + ], + [ + -0.2339, + -1.0397 + ], + [ + 0.0442, + -1.028 + ], + [ + 0.2982, + -1.0291 + ], + [ + 0.553, + -1.067 + ], + [ + 0.8784, + -1.0285 + ], + [ + 1.125, + -1.0434 + ], + [ + 1.4259, + -1.0245 + ], + [ + 1.7058, + -1.05 + ], + [ + 1.973, + -1.0492 + ], + [ + 2.2685, + -1.0245 + ], + [ + -2.3626, + -0.8167 + ], + [ + -2.0947, + -0.7888 + ], + [ + -1.7823, + -0.8247 + ], + [ + -1.5443, + -0.7973 + ], + [ + -1.2598, + -0.7844 + ], + [ + -0.9706, + -0.8225 + ], + [ + -0.6564, + -0.8213 + ], + [ + -0.4121, + -0.822 + ], + [ + -0.1415, + -0.8294 + ], + [ + 0.1351, + -0.8208 + ], + [ + 0.415, + -0.7975 + ], + [ + 0.7301, + -0.8085 + ], + [ + 0.9898, + -0.8023 + ], + [ + 1.301, + -0.7829 + ], + [ + 1.5344, + -0.8204 + ], + [ + 1.8253, + -0.8166 + ], + [ + 2.124, + -0.8073 + ], + [ + 2.3748, + -0.8097 + ], + [ + -2.4987, + -0.5824 + ], + [ + -2.2013, + -0.5852 + ], + [ + -1.9162, + -0.5826 + ], + [ + -1.6408, + -0.5376 + ], + [ + -1.3531, + -0.5432 + ], + [ + -1.0843, + -0.5521 + ], + [ + -0.8044, + -0.5806 + ], + [ + -0.538, + -0.5593 + ], + [ + -0.24, + -0.5622 + ], + [ + 0.0136, + -0.5523 + ], + [ + 0.2869, + -0.5803 + ], + [ + 0.5788, + -0.5648 + ], + [ + 0.845, + -0.5675 + ], + [ + 1.1325, + -0.5698 + ], + [ + 1.4132, + -0.5497 + ], + [ + 1.6886, + -0.535 + ], + [ + 2.0033, + -0.5612 + ], + [ + 2.2504, + -0.5581 + ], + [ + -2.3405, + -0.3091 + ], + [ + -2.063, + -0.3158 + ], + [ + -1.7926, + -0.3295 + ], + [ + -1.5068, + -0.3343 + ], + [ + -1.2291, + -0.2974 + ], + [ + -0.9806, + -0.3112 + ], + [ + -0.7027, + -0.305 + ], + [ + -0.4233, + -0.2932 + ], + [ + -0.1403, + -0.2919 + ], + [ + 0.1768, + -0.3124 + ], + [ + 0.4558, + -0.3011 + ], + [ + 0.745, + -0.3314 + ], + [ + 1.005, + -0.3403 + ], + [ + 1.2865, + -0.336 + ], + [ + 1.5636, + -0.3136 + ], + [ + 1.8381, + -0.3163 + ], + [ + 2.1348, + -0.3009 + ], + [ + 2.3996, + -0.312 + ], + [ + -2.4759, + -0.0964 + ], + [ + -2.2414, + -0.0982 + ], + [ + -1.9312, + -0.0797 + ], + [ + -1.6446, + -0.0655 + ], + [ + -1.3893, + -0.0528 + ], + [ + -1.0853, + -0.088 + ], + [ + -0.8276, + -0.0855 + ], + [ + -0.5592, + -0.0949 + ], + [ + -0.2356, + -0.0786 + ], + [ + 0.0135, + -0.0623 + ], + [ + 0.303, + -0.0507 + ], + [ + 0.5957, + -0.0763 + ], + [ + 0.8531, + -0.0479 + ], + [ + 1.1522, + -0.0499 + ], + [ + 1.3986, + -0.0555 + ], + [ + 1.7077, + -0.0963 + ], + [ + 1.9849, + -0.0647 + ], + [ + 2.2327, + -0.0777 + ], + [ + -2.3418, + 0.1559 + ], + [ + -2.0823, + 0.1641 + ], + [ + -1.8111, + 0.1908 + ], + [ + -1.5042, + 0.1456 + ], + [ + -1.2122, + 0.1886 + ], + [ + -0.9721, + 0.1862 + ], + [ + -0.702, + 0.1953 + ], + [ + -0.3907, + 0.1758 + ], + [ + -0.1429, + 0.1896 + ], + [ + 0.1336, + 0.1528 + ], + [ + 0.42, + 0.1829 + ], + [ + 0.7031, + 0.1904 + ], + [ + 1.0088, + 0.1414 + ], + [ + 1.2523, + 0.1423 + ], + [ + 1.5659, + 0.1842 + ], + [ + 1.8254, + 0.1869 + ], + [ + 2.0952, + 0.1842 + ], + [ + 2.424, + 0.1826 + ], + [ + -2.4889, + 0.4288 + ], + [ + -2.2458, + 0.3931 + ], + [ + -1.961, + 0.4101 + ], + [ + -1.6463, + 0.4171 + ], + [ + -1.3603, + 0.3905 + ], + [ + -1.0869, + 0.3927 + ], + [ + -0.8328, + 0.4216 + ], + [ + -0.5131, + 0.4161 + ], + [ + -0.2849, + 0.4164 + ], + [ + -0.0056, + 0.4314 + ], + [ + 0.3117, + 0.3915 + ], + [ + 0.5571, + 0.3921 + ], + [ + 0.8869, + 0.4075 + ], + [ + 1.1559, + 0.4175 + ], + [ + 1.4241, + 0.39 + ], + [ + 1.725, + 0.3987 + ], + [ + 1.9844, + 0.421 + ], + [ + 2.2684, + 0.4345 + ], + [ + -2.3797, + 0.6528 + ], + [ + -2.0854, + 0.6509 + ], + [ + -1.8213, + 0.6318 + ], + [ + -1.5324, + 0.6414 + ], + [ + -1.244, + 0.6585 + ], + [ + -0.9525, + 0.6474 + ], + [ + -0.6851, + 0.6365 + ], + [ + -0.3951, + 0.6421 + ], + [ + -0.146, + 0.6478 + ], + [ + 0.1586, + 0.637 + ], + [ + 0.4441, + 0.656 + ], + [ + 0.7313, + 0.6606 + ], + [ + 1.0085, + 0.642 + ], + [ + 1.2961, + 0.6551 + ], + [ + 1.5562, + 0.6594 + ], + [ + 1.8322, + 0.653 + ], + [ + 2.1333, + 0.674 + ], + [ + 2.4236, + 0.6525 + ], + [ + -2.4989, + 0.9116 + ], + [ + -2.2304, + 0.9137 + ], + [ + -1.9403, + 0.8733 + ], + [ + -1.684, + 0.914 + ], + [ + -1.4049, + 0.8825 + ], + [ + -1.1093, + 0.8765 + ], + [ + -0.8304, + 0.9084 + ], + [ + -0.5672, + 0.9131 + ], + [ + -0.24, + 0.8877 + ], + [ + 0.0006, + 0.9005 + ], + [ + 0.2787, + 0.8873 + ], + [ + 0.6057, + 0.9226 + ], + [ + 0.8752, + 0.8842 + ], + [ + 1.1505, + 0.9063 + ], + [ + 1.4137, + 0.9027 + ], + [ + 1.6726, + 0.8785 + ], + [ + 1.9814, + 0.876 + ], + [ + 2.2413, + 0.9137 + ], + [ + -2.3326, + 1.1404 + ], + [ + -2.061, + 1.1648 + ], + [ + -1.8201, + 1.1344 + ], + [ + -1.526, + 1.1138 + ], + [ + -1.2257, + 1.1336 + ], + [ + -0.9617, + 1.1177 + ], + [ + -0.6979, + 1.1601 + ], + [ + -0.4255, + 1.1223 + ], + [ + -0.1316, + 1.1368 + ], + [ + 0.1648, + 1.1369 + ], + [ + 0.4167, + 1.123 + ], + [ + 0.7392, + 1.145 + ], + [ + 1.0084, + 1.1468 + ], + [ + 1.2947, + 1.1126 + ], + [ + 1.5525, + 1.1395 + ], + [ + 1.831, + 1.1566 + ], + [ + 2.119, + 1.1523 + ], + [ + 2.4197, + 1.1376 + ], + [ + -2.4771, + 1.3847 + ], + [ + -2.2004, + 1.3709 + ], + [ + -1.9401, + 1.3816 + ], + [ + -1.6821, + 1.3741 + ], + [ + -1.3566, + 1.3871 + ], + [ + -1.1181, + 1.3708 + ], + [ + -0.8373, + 1.3532 + ], + [ + -0.5161, + 1.3769 + ], + [ + -0.2708, + 1.3853 + ], + [ + -0.0076, + 1.3674 + ], + [ + 0.3114, + 1.3873 + ], + [ + 0.607, + 1.3865 + ], + [ + 0.8587, + 1.3944 + ], + [ + 1.1626, + 1.3922 + ], + [ + 1.4459, + 1.3956 + ], + [ + 1.7205, + 1.3582 + ], + [ + 1.993, + 1.3764 + ], + [ + 2.263, + 1.3884 + ], + [ + -2.3337, + 1.6494 + ], + [ + -2.0919, + 1.6354 + ], + [ + -1.786, + 1.6137 + ], + [ + -1.5411, + 1.5966 + ], + [ + -1.2245, + 1.6217 + ], + [ + -0.9328, + 1.6203 + ], + [ + -0.6532, + 1.6173 + ], + [ + -0.3836, + 1.599 + ], + [ + -0.1169, + 1.6392 + ], + [ + 0.1838, + 1.6403 + ], + [ + 0.4141, + 1.6152 + ], + [ + 0.6947, + 1.6004 + ], + [ + 1.0098, + 1.6342 + ], + [ + 1.2953, + 1.6427 + ], + [ + 1.5734, + 1.6391 + ], + [ + 1.8147, + 1.6074 + ], + [ + 2.1268, + 1.6423 + ], + [ + 2.3723, + 1.6231 + ], + [ + -2.4901, + 1.8384 + ], + [ + -2.2255, + 1.8869 + ], + [ + -1.9304, + 1.8501 + ], + [ + -1.6402, + 1.8563 + ], + [ + -1.3602, + 1.8535 + ], + [ + -1.0949, + 1.859 + ], + [ + -0.8326, + 1.8864 + ], + [ + -0.5575, + 1.8415 + ], + [ + -0.2689, + 1.877 + ], + [ + 0.0372, + 1.8927 + ], + [ + 0.2886, + 1.8596 + ], + [ + 0.5597, + 1.869 + ], + [ + 0.8879, + 1.876 + ], + [ + 1.1453, + 1.8587 + ], + [ + 1.4433, + 1.8646 + ], + [ + 1.6795, + 1.8572 + ], + [ + 1.9558, + 1.8481 + ], + [ + 2.233, + 1.8622 + ], + [ + -2.3525, + 2.0985 + ], + [ + -2.0845, + 2.133 + ], + [ + -1.7859, + 2.1095 + ], + [ + -1.5321, + 2.1295 + ], + [ + -1.2548, + 2.0975 + ], + [ + -0.9371, + 2.1089 + ], + [ + -0.6664, + 2.1123 + ], + [ + -0.3914, + 2.096 + ], + [ + -0.1345, + 2.0973 + ], + [ + 0.1407, + 2.1282 + ], + [ + 0.4279, + 2.1107 + ], + [ + 0.7364, + 2.1231 + ], + [ + 0.9965, + 2.1059 + ], + [ + 1.3077, + 2.117 + ], + [ + 1.5776, + 2.1298 + ], + [ + 1.8561, + 2.0896 + ], + [ + 2.1235, + 2.085 + ], + [ + 2.4086, + 2.1328 + ], + [ + -2.4993, + 2.346 + ], + [ + -2.246, + 2.3755 + ], + [ + -1.9622, + 2.324 + ], + [ + -1.6742, + 2.3254 + ], + [ + -1.3825, + 2.3506 + ], + [ + -1.1105, + 2.3246 + ], + [ + -0.8418, + 2.3433 + ], + [ + -0.5646, + 2.3608 + ], + [ + -0.2764, + 2.3387 + ], + [ + 0.0141, + 2.3451 + ], + [ + 0.2721, + 2.328 + ], + [ + 0.6003, + 2.3218 + ], + [ + 0.8605, + 2.3492 + ], + [ + 1.1307, + 2.3459 + ], + [ + 1.4357, + 2.3688 + ], + [ + 1.6866, + 2.3398 + ], + [ + 1.9656, + 2.3486 + ], + [ + 2.2703, + 2.3345 + ], + [ + -2.3695, + 2.6163 + ], + [ + -2.1053, + 2.59 + ], + [ + -1.7882, + 2.5727 + ], + [ + -1.5453, + 2.572 + ], + [ + -1.2165, + 2.5647 + ], + [ + -0.9775, + 2.566 + ], + [ + -0.7018, + 2.599 + ], + [ + -0.4145, + 2.5961 + ], + [ + -0.1149, + 2.6118 + ], + [ + 0.1323, + 2.612 + ], + [ + 0.4467, + 2.5733 + ], + [ + 0.7353, + 2.6121 + ], + [ + 0.9862, + 2.6157 + ], + [ + 1.2774, + 2.598 + ], + [ + 1.5871, + 2.5845 + ], + [ + 1.8575, + 2.5821 + ], + [ + 2.1368, + 2.5979 + ], + [ + 2.3841, + 2.5874 + ] + ], + "triangles": [ + [ + 107, + 89, + 71 + ], + [ + 252, + 360, + 108 + ], + [ + 0, + 72, + 108 + ], + [ + 216, + 252, + 108 + ], + [ + 107, + 179, + 143 + ], + [ + 380, + 363, + 381 + ], + [ + 251, + 250, + 233 + ], + [ + 370, + 387, + 369 + ], + [ + 36, + 37, + 54 + ], + [ + 36, + 72, + 0 + ], + [ + 72, + 36, + 54 + ], + [ + 37, + 55, + 54 + ], + [ + 18, + 36, + 0 + ], + [ + 36, + 18, + 37 + ], + [ + 198, + 216, + 180 + ], + [ + 216, + 198, + 217 + ], + [ + 234, + 216, + 217 + ], + [ + 234, + 253, + 252 + ], + [ + 216, + 234, + 252 + ], + [ + 150, + 149, + 131 + ], + [ + 20, + 21, + 39 + ], + [ + 21, + 20, + 3 + ], + [ + 21, + 40, + 39 + ], + [ + 40, + 21, + 22 + ], + [ + 40, + 41, + 58 + ], + [ + 41, + 40, + 22 + ], + [ + 40, + 57, + 39 + ], + [ + 57, + 40, + 58 + ], + [ + 144, + 216, + 108 + ], + [ + 216, + 144, + 180 + ], + [ + 146, + 163, + 145 + ], + [ + 364, + 363, + 345 + ], + [ + 363, + 364, + 381 + ], + [ + 381, + 384, + 378 + ], + [ + 215, + 179, + 107 + ], + [ + 215, + 107, + 71 + ], + [ + 215, + 251, + 233 + ], + [ + 161, + 179, + 178 + ], + [ + 160, + 161, + 178 + ], + [ + 179, + 161, + 143 + ], + [ + 161, + 142, + 143 + ], + [ + 142, + 161, + 160 + ], + [ + 177, + 160, + 178 + ], + [ + 270, + 289, + 288 + ], + [ + 253, + 270, + 252 + ], + [ + 270, + 288, + 252 + ], + [ + 307, + 289, + 290 + ], + [ + 308, + 307, + 290 + ], + [ + 307, + 308, + 326 + ], + [ + 308, + 327, + 326 + ], + [ + 327, + 308, + 309 + ], + [ + 362, + 363, + 380 + ], + [ + 212, + 194, + 195 + ], + [ + 175, + 194, + 193 + ], + [ + 250, + 232, + 233 + ], + [ + 267, + 268, + 285 + ], + [ + 250, + 269, + 268 + ], + [ + 269, + 250, + 251 + ], + [ + 287, + 269, + 251 + ], + [ + 305, + 287, + 323 + ], + [ + 301, + 302, + 319 + ], + [ + 287, + 359, + 323 + ], + [ + 370, + 352, + 371 + ], + [ + 394, + 390, + 391 + ], + [ + 391, + 390, + 373 + ], + [ + 389, + 390, + 387 + ], + [ + 388, + 370, + 371 + ], + [ + 389, + 388, + 371 + ], + [ + 370, + 388, + 387 + ], + [ + 388, + 389, + 387 + ], + [ + 372, + 389, + 371 + ], + [ + 390, + 372, + 373 + ], + [ + 372, + 390, + 389 + ], + [ + 38, + 55, + 37 + ], + [ + 38, + 20, + 39 + ], + [ + 20, + 2, + 3 + ], + [ + 35, + 34, + 17 + ], + [ + 181, + 198, + 180 + ], + [ + 181, + 163, + 182 + ], + [ + 235, + 234, + 217 + ], + [ + 254, + 235, + 236 + ], + [ + 235, + 254, + 253 + ], + [ + 234, + 235, + 253 + ], + [ + 149, + 130, + 131 + ], + [ + 168, + 167, + 150 + ], + [ + 167, + 149, + 150 + ], + [ + 166, + 167, + 185 + ], + [ + 167, + 166, + 149 + ], + [ + 227, + 226, + 208 + ], + [ + 225, + 226, + 243 + ], + [ + 245, + 227, + 228 + ], + [ + 5, + 23, + 22 + ], + [ + 23, + 41, + 22 + ], + [ + 23, + 5, + 6 + ], + [ + 4, + 5, + 22 + ], + [ + 21, + 4, + 22 + ], + [ + 4, + 21, + 3 + ], + [ + 4, + 7, + 6 + ], + [ + 5, + 4, + 6 + ], + [ + 8, + 4, + 0 + ], + [ + 7, + 4, + 8 + ], + [ + 76, + 57, + 58 + ], + [ + 126, + 144, + 108 + ], + [ + 144, + 126, + 145 + ], + [ + 72, + 90, + 108 + ], + [ + 55, + 73, + 54 + ], + [ + 73, + 72, + 54 + ], + [ + 73, + 90, + 72 + ], + [ + 90, + 73, + 91 + ], + [ + 164, + 163, + 146 + ], + [ + 163, + 164, + 182 + ], + [ + 255, + 254, + 236 + ], + [ + 255, + 256, + 273 + ], + [ + 272, + 255, + 273 + ], + [ + 255, + 272, + 254 + ], + [ + 309, + 292, + 310 + ], + [ + 292, + 293, + 310 + ], + [ + 240, + 222, + 223 + ], + [ + 187, + 168, + 169 + ], + [ + 188, + 187, + 169 + ], + [ + 206, + 188, + 189 + ], + [ + 242, + 225, + 243 + ], + [ + 364, + 382, + 381 + ], + [ + 382, + 364, + 365 + ], + [ + 382, + 384, + 381 + ], + [ + 346, + 364, + 345 + ], + [ + 364, + 346, + 365 + ], + [ + 346, + 347, + 365 + ], + [ + 347, + 346, + 329 + ], + [ + 125, + 107, + 143 + ], + [ + 142, + 125, + 143 + ], + [ + 177, + 159, + 160 + ], + [ + 291, + 272, + 273 + ], + [ + 292, + 291, + 273 + ], + [ + 291, + 292, + 309 + ], + [ + 272, + 291, + 290 + ], + [ + 291, + 308, + 290 + ], + [ + 308, + 291, + 309 + ], + [ + 270, + 271, + 289 + ], + [ + 289, + 271, + 290 + ], + [ + 271, + 272, + 290 + ], + [ + 271, + 270, + 253 + ], + [ + 254, + 271, + 253 + ], + [ + 272, + 271, + 254 + ], + [ + 325, + 307, + 326 + ], + [ + 289, + 306, + 288 + ], + [ + 307, + 306, + 289 + ], + [ + 325, + 306, + 307 + ], + [ + 362, + 344, + 363 + ], + [ + 363, + 344, + 345 + ], + [ + 344, + 327, + 345 + ], + [ + 327, + 344, + 326 + ], + [ + 379, + 362, + 380 + ], + [ + 379, + 381, + 378 + ], + [ + 379, + 380, + 381 + ], + [ + 325, + 343, + 342 + ], + [ + 343, + 325, + 326 + ], + [ + 344, + 343, + 326 + ], + [ + 343, + 344, + 362 + ], + [ + 227, + 209, + 228 + ], + [ + 192, + 209, + 191 + ], + [ + 209, + 208, + 191 + ], + [ + 209, + 227, + 208 + ], + [ + 210, + 229, + 228 + ], + [ + 210, + 192, + 193 + ], + [ + 209, + 210, + 228 + ], + [ + 210, + 209, + 192 + ], + [ + 248, + 267, + 266 + ], + [ + 194, + 176, + 195 + ], + [ + 176, + 194, + 175 + ], + [ + 176, + 177, + 195 + ], + [ + 176, + 159, + 177 + ], + [ + 176, + 175, + 158 + ], + [ + 159, + 176, + 158 + ], + [ + 211, + 194, + 212 + ], + [ + 194, + 211, + 193 + ], + [ + 211, + 210, + 193 + ], + [ + 210, + 211, + 229 + ], + [ + 213, + 212, + 195 + ], + [ + 214, + 215, + 233 + ], + [ + 232, + 214, + 233 + ], + [ + 213, + 214, + 232 + ], + [ + 284, + 267, + 285 + ], + [ + 267, + 284, + 266 + ], + [ + 302, + 320, + 319 + ], + [ + 337, + 318, + 319 + ], + [ + 318, + 301, + 319 + ], + [ + 359, + 341, + 323 + ], + [ + 394, + 377, + 395 + ], + [ + 377, + 359, + 395 + ], + [ + 392, + 394, + 391 + ], + [ + 334, + 333, + 315 + ], + [ + 336, + 318, + 337 + ], + [ + 318, + 336, + 317 + ], + [ + 38, + 56, + 55 + ], + [ + 57, + 56, + 39 + ], + [ + 56, + 38, + 39 + ], + [ + 1, + 18, + 0 + ], + [ + 2, + 1, + 3 + ], + [ + 4, + 1, + 0 + ], + [ + 1, + 4, + 3 + ], + [ + 38, + 19, + 20 + ], + [ + 19, + 2, + 20 + ], + [ + 19, + 38, + 37 + ], + [ + 19, + 1, + 2 + ], + [ + 18, + 19, + 37 + ], + [ + 1, + 19, + 18 + ], + [ + 53, + 35, + 71 + ], + [ + 53, + 34, + 35 + ], + [ + 16, + 13, + 17 + ], + [ + 34, + 16, + 17 + ], + [ + 14, + 15, + 32 + ], + [ + 16, + 14, + 13 + ], + [ + 14, + 16, + 15 + ], + [ + 15, + 33, + 32 + ], + [ + 33, + 16, + 34 + ], + [ + 16, + 33, + 15 + ], + [ + 163, + 162, + 145 + ], + [ + 181, + 162, + 163 + ], + [ + 162, + 144, + 145 + ], + [ + 144, + 162, + 180 + ], + [ + 162, + 181, + 180 + ], + [ + 218, + 235, + 217 + ], + [ + 235, + 218, + 236 + ], + [ + 218, + 219, + 236 + ], + [ + 219, + 218, + 200 + ], + [ + 148, + 130, + 149 + ], + [ + 148, + 166, + 165 + ], + [ + 166, + 148, + 149 + ], + [ + 229, + 246, + 228 + ], + [ + 246, + 245, + 228 + ], + [ + 245, + 244, + 227 + ], + [ + 244, + 262, + 243 + ], + [ + 226, + 244, + 243 + ], + [ + 244, + 226, + 227 + ], + [ + 175, + 157, + 158 + ], + [ + 25, + 7, + 8 + ], + [ + 168, + 151, + 169 + ], + [ + 151, + 168, + 150 + ], + [ + 130, + 113, + 131 + ], + [ + 118, + 99, + 100 + ], + [ + 25, + 44, + 43 + ], + [ + 64, + 46, + 47 + ], + [ + 46, + 64, + 63 + ], + [ + 23, + 42, + 41 + ], + [ + 60, + 42, + 43 + ], + [ + 75, + 76, + 93 + ], + [ + 76, + 75, + 57 + ], + [ + 75, + 56, + 57 + ], + [ + 92, + 110, + 91 + ], + [ + 92, + 75, + 93 + ], + [ + 148, + 129, + 130 + ], + [ + 109, + 90, + 91 + ], + [ + 110, + 109, + 91 + ], + [ + 109, + 126, + 108 + ], + [ + 90, + 109, + 108 + ], + [ + 127, + 146, + 145 + ], + [ + 127, + 109, + 110 + ], + [ + 126, + 127, + 145 + ], + [ + 109, + 127, + 126 + ], + [ + 256, + 237, + 238 + ], + [ + 219, + 237, + 236 + ], + [ + 237, + 255, + 236 + ], + [ + 255, + 237, + 256 + ], + [ + 201, + 219, + 200 + ], + [ + 274, + 275, + 293 + ], + [ + 274, + 292, + 273 + ], + [ + 292, + 274, + 293 + ], + [ + 256, + 274, + 273 + ], + [ + 240, + 239, + 222 + ], + [ + 186, + 203, + 185 + ], + [ + 187, + 186, + 168 + ], + [ + 167, + 186, + 185 + ], + [ + 186, + 167, + 168 + ], + [ + 224, + 206, + 225 + ], + [ + 242, + 224, + 225 + ], + [ + 205, + 187, + 188 + ], + [ + 206, + 205, + 188 + ], + [ + 205, + 224, + 223 + ], + [ + 224, + 205, + 206 + ], + [ + 207, + 226, + 225 + ], + [ + 206, + 207, + 225 + ], + [ + 207, + 206, + 189 + ], + [ + 226, + 207, + 208 + ], + [ + 296, + 314, + 313 + ], + [ + 333, + 314, + 315 + ], + [ + 314, + 332, + 313 + ], + [ + 314, + 333, + 332 + ], + [ + 241, + 242, + 260 + ], + [ + 259, + 241, + 260 + ], + [ + 241, + 259, + 240 + ], + [ + 241, + 240, + 223 + ], + [ + 224, + 241, + 223 + ], + [ + 241, + 224, + 242 + ], + [ + 383, + 382, + 365 + ], + [ + 382, + 383, + 384 + ], + [ + 332, + 331, + 313 + ], + [ + 331, + 312, + 313 + ], + [ + 294, + 275, + 276 + ], + [ + 275, + 294, + 293 + ], + [ + 295, + 296, + 313 + ], + [ + 312, + 295, + 313 + ], + [ + 295, + 294, + 276 + ], + [ + 294, + 295, + 312 + ], + [ + 327, + 328, + 345 + ], + [ + 328, + 346, + 345 + ], + [ + 328, + 309, + 310 + ], + [ + 328, + 327, + 309 + ], + [ + 329, + 328, + 310 + ], + [ + 346, + 328, + 329 + ], + [ + 331, + 349, + 348 + ], + [ + 349, + 331, + 332 + ], + [ + 349, + 367, + 348 + ], + [ + 367, + 349, + 368 + ], + [ + 384, + 386, + 378 + ], + [ + 387, + 386, + 369 + ], + [ + 386, + 368, + 369 + ], + [ + 386, + 390, + 378 + ], + [ + 390, + 386, + 387 + ], + [ + 351, + 352, + 370 + ], + [ + 351, + 370, + 369 + ], + [ + 351, + 334, + 352 + ], + [ + 334, + 351, + 333 + ], + [ + 324, + 325, + 342 + ], + [ + 324, + 306, + 325 + ], + [ + 306, + 324, + 288 + ], + [ + 360, + 324, + 342 + ], + [ + 324, + 360, + 252 + ], + [ + 288, + 324, + 252 + ], + [ + 379, + 361, + 362 + ], + [ + 361, + 343, + 362 + ], + [ + 360, + 361, + 378 + ], + [ + 361, + 379, + 378 + ], + [ + 361, + 360, + 342 + ], + [ + 343, + 361, + 342 + ], + [ + 230, + 211, + 212 + ], + [ + 211, + 230, + 229 + ], + [ + 231, + 213, + 232 + ], + [ + 230, + 231, + 248 + ], + [ + 213, + 231, + 212 + ], + [ + 231, + 230, + 212 + ], + [ + 214, + 197, + 215 + ], + [ + 179, + 197, + 178 + ], + [ + 215, + 197, + 179 + ], + [ + 196, + 213, + 195 + ], + [ + 196, + 214, + 213 + ], + [ + 196, + 197, + 214 + ], + [ + 177, + 196, + 195 + ], + [ + 196, + 177, + 178 + ], + [ + 197, + 196, + 178 + ], + [ + 303, + 284, + 285 + ], + [ + 284, + 303, + 302 + ], + [ + 303, + 320, + 302 + ], + [ + 304, + 303, + 285 + ], + [ + 284, + 283, + 266 + ], + [ + 283, + 301, + 282 + ], + [ + 301, + 283, + 302 + ], + [ + 283, + 284, + 302 + ], + [ + 304, + 286, + 305 + ], + [ + 286, + 269, + 287 + ], + [ + 305, + 286, + 287 + ], + [ + 286, + 304, + 285 + ], + [ + 268, + 286, + 285 + ], + [ + 269, + 286, + 268 + ], + [ + 322, + 305, + 323 + ], + [ + 322, + 304, + 305 + ], + [ + 341, + 322, + 323 + ], + [ + 322, + 341, + 340 + ], + [ + 355, + 338, + 356 + ], + [ + 320, + 338, + 319 + ], + [ + 338, + 337, + 319 + ], + [ + 338, + 355, + 337 + ], + [ + 300, + 318, + 317 + ], + [ + 300, + 281, + 282 + ], + [ + 301, + 300, + 282 + ], + [ + 318, + 300, + 301 + ], + [ + 281, + 264, + 282 + ], + [ + 246, + 264, + 245 + ], + [ + 263, + 281, + 280 + ], + [ + 263, + 244, + 245 + ], + [ + 264, + 263, + 245 + ], + [ + 263, + 264, + 281 + ], + [ + 262, + 263, + 280 + ], + [ + 244, + 263, + 262 + ], + [ + 299, + 300, + 317 + ], + [ + 300, + 299, + 281 + ], + [ + 299, + 298, + 280 + ], + [ + 281, + 299, + 280 + ], + [ + 392, + 393, + 394 + ], + [ + 358, + 341, + 359 + ], + [ + 377, + 358, + 359 + ], + [ + 341, + 358, + 340 + ], + [ + 358, + 357, + 340 + ], + [ + 374, + 355, + 356 + ], + [ + 374, + 392, + 391 + ], + [ + 374, + 391, + 373 + ], + [ + 355, + 374, + 373 + ], + [ + 316, + 334, + 315 + ], + [ + 298, + 316, + 315 + ], + [ + 316, + 299, + 317 + ], + [ + 299, + 316, + 298 + ], + [ + 336, + 335, + 317 + ], + [ + 334, + 335, + 352 + ], + [ + 335, + 316, + 317 + ], + [ + 316, + 335, + 334 + ], + [ + 372, + 354, + 373 + ], + [ + 354, + 336, + 337 + ], + [ + 354, + 355, + 373 + ], + [ + 355, + 354, + 337 + ], + [ + 73, + 74, + 91 + ], + [ + 74, + 92, + 91 + ], + [ + 74, + 73, + 55 + ], + [ + 56, + 74, + 55 + ], + [ + 75, + 74, + 56 + ], + [ + 92, + 74, + 75 + ], + [ + 89, + 70, + 71 + ], + [ + 70, + 53, + 71 + ], + [ + 31, + 14, + 32 + ], + [ + 14, + 31, + 13 + ], + [ + 198, + 199, + 217 + ], + [ + 199, + 218, + 217 + ], + [ + 218, + 199, + 200 + ], + [ + 200, + 199, + 182 + ], + [ + 199, + 181, + 182 + ], + [ + 181, + 199, + 198 + ], + [ + 247, + 246, + 229 + ], + [ + 230, + 247, + 229 + ], + [ + 247, + 230, + 248 + ], + [ + 247, + 248, + 266 + ], + [ + 208, + 190, + 191 + ], + [ + 207, + 190, + 208 + ], + [ + 171, + 190, + 189 + ], + [ + 190, + 207, + 189 + ], + [ + 119, + 118, + 100 + ], + [ + 120, + 121, + 138 + ], + [ + 122, + 121, + 103 + ], + [ + 121, + 102, + 103 + ], + [ + 102, + 121, + 120 + ], + [ + 141, + 142, + 160 + ], + [ + 159, + 141, + 160 + ], + [ + 172, + 171, + 154 + ], + [ + 155, + 172, + 154 + ], + [ + 190, + 172, + 191 + ], + [ + 172, + 190, + 171 + ], + [ + 137, + 120, + 138 + ], + [ + 137, + 119, + 120 + ], + [ + 173, + 192, + 191 + ], + [ + 172, + 173, + 191 + ], + [ + 173, + 172, + 155 + ], + [ + 7, + 24, + 6 + ], + [ + 25, + 24, + 7 + ], + [ + 24, + 23, + 6 + ], + [ + 24, + 42, + 23 + ], + [ + 24, + 25, + 43 + ], + [ + 42, + 24, + 43 + ], + [ + 117, + 99, + 118 + ], + [ + 99, + 117, + 98 + ], + [ + 171, + 153, + 154 + ], + [ + 113, + 114, + 131 + ], + [ + 114, + 113, + 95 + ], + [ + 113, + 94, + 95 + ], + [ + 76, + 94, + 93 + ], + [ + 29, + 48, + 47 + ], + [ + 64, + 82, + 63 + ], + [ + 99, + 82, + 100 + ], + [ + 26, + 25, + 8 + ], + [ + 26, + 44, + 25 + ], + [ + 11, + 28, + 10 + ], + [ + 28, + 11, + 29 + ], + [ + 28, + 29, + 47 + ], + [ + 46, + 28, + 47 + ], + [ + 59, + 42, + 60 + ], + [ + 78, + 59, + 60 + ], + [ + 41, + 59, + 58 + ], + [ + 42, + 59, + 41 + ], + [ + 127, + 128, + 146 + ], + [ + 128, + 127, + 110 + ], + [ + 164, + 147, + 165 + ], + [ + 147, + 164, + 146 + ], + [ + 147, + 148, + 165 + ], + [ + 147, + 129, + 148 + ], + [ + 128, + 147, + 146 + ], + [ + 147, + 128, + 129 + ], + [ + 94, + 112, + 93 + ], + [ + 112, + 94, + 113 + ], + [ + 112, + 113, + 130 + ], + [ + 129, + 112, + 130 + ], + [ + 183, + 200, + 182 + ], + [ + 183, + 201, + 200 + ], + [ + 164, + 183, + 182 + ], + [ + 183, + 164, + 165 + ], + [ + 201, + 220, + 219 + ], + [ + 237, + 220, + 238 + ], + [ + 220, + 237, + 219 + ], + [ + 274, + 257, + 275 + ], + [ + 239, + 257, + 238 + ], + [ + 257, + 256, + 238 + ], + [ + 257, + 274, + 256 + ], + [ + 259, + 258, + 240 + ], + [ + 258, + 239, + 240 + ], + [ + 275, + 258, + 276 + ], + [ + 258, + 259, + 276 + ], + [ + 257, + 258, + 275 + ], + [ + 258, + 257, + 239 + ], + [ + 186, + 204, + 203 + ], + [ + 204, + 186, + 187 + ], + [ + 205, + 204, + 187 + ], + [ + 203, + 204, + 222 + ], + [ + 222, + 204, + 223 + ], + [ + 204, + 205, + 223 + ], + [ + 261, + 278, + 260 + ], + [ + 262, + 261, + 243 + ], + [ + 261, + 242, + 243 + ], + [ + 242, + 261, + 260 + ], + [ + 297, + 298, + 315 + ], + [ + 314, + 297, + 315 + ], + [ + 297, + 314, + 296 + ], + [ + 278, + 297, + 296 + ], + [ + 347, + 366, + 365 + ], + [ + 366, + 383, + 365 + ], + [ + 366, + 347, + 348 + ], + [ + 367, + 366, + 348 + ], + [ + 383, + 366, + 384 + ], + [ + 366, + 367, + 384 + ], + [ + 347, + 330, + 348 + ], + [ + 330, + 331, + 348 + ], + [ + 330, + 347, + 329 + ], + [ + 331, + 330, + 312 + ], + [ + 277, + 259, + 260 + ], + [ + 278, + 277, + 260 + ], + [ + 259, + 277, + 276 + ], + [ + 277, + 295, + 276 + ], + [ + 277, + 278, + 296 + ], + [ + 295, + 277, + 296 + ], + [ + 333, + 350, + 332 + ], + [ + 350, + 349, + 332 + ], + [ + 351, + 350, + 333 + ], + [ + 350, + 351, + 369 + ], + [ + 368, + 350, + 369 + ], + [ + 349, + 350, + 368 + ], + [ + 385, + 367, + 368 + ], + [ + 386, + 385, + 368 + ], + [ + 367, + 385, + 384 + ], + [ + 385, + 386, + 384 + ], + [ + 249, + 231, + 232 + ], + [ + 249, + 250, + 268 + ], + [ + 249, + 232, + 250 + ], + [ + 231, + 249, + 248 + ], + [ + 267, + 249, + 268 + ], + [ + 248, + 249, + 267 + ], + [ + 321, + 303, + 304 + ], + [ + 303, + 321, + 320 + ], + [ + 321, + 322, + 340 + ], + [ + 322, + 321, + 304 + ], + [ + 283, + 265, + 266 + ], + [ + 265, + 247, + 266 + ], + [ + 265, + 283, + 282 + ], + [ + 264, + 265, + 282 + ], + [ + 265, + 264, + 246 + ], + [ + 247, + 265, + 246 + ], + [ + 376, + 377, + 394 + ], + [ + 393, + 376, + 394 + ], + [ + 376, + 358, + 377 + ], + [ + 358, + 376, + 357 + ], + [ + 338, + 339, + 356 + ], + [ + 339, + 357, + 356 + ], + [ + 357, + 339, + 340 + ], + [ + 339, + 338, + 320 + ], + [ + 339, + 321, + 340 + ], + [ + 321, + 339, + 320 + ], + [ + 375, + 374, + 356 + ], + [ + 374, + 375, + 392 + ], + [ + 357, + 375, + 356 + ], + [ + 376, + 375, + 357 + ], + [ + 375, + 393, + 392 + ], + [ + 375, + 376, + 393 + ], + [ + 353, + 335, + 336 + ], + [ + 354, + 353, + 336 + ], + [ + 352, + 353, + 371 + ], + [ + 335, + 353, + 352 + ], + [ + 353, + 372, + 371 + ], + [ + 353, + 354, + 372 + ], + [ + 52, + 33, + 34 + ], + [ + 53, + 52, + 34 + ], + [ + 70, + 52, + 53 + ], + [ + 69, + 52, + 70 + ], + [ + 33, + 51, + 32 + ], + [ + 52, + 51, + 33 + ], + [ + 51, + 52, + 69 + ], + [ + 139, + 121, + 122 + ], + [ + 157, + 139, + 158 + ], + [ + 139, + 157, + 138 + ], + [ + 121, + 139, + 138 + ], + [ + 124, + 125, + 142 + ], + [ + 141, + 124, + 142 + ], + [ + 119, + 136, + 118 + ], + [ + 137, + 136, + 119 + ], + [ + 136, + 137, + 155 + ], + [ + 136, + 155, + 154 + ], + [ + 192, + 174, + 193 + ], + [ + 173, + 174, + 192 + ], + [ + 174, + 175, + 193 + ], + [ + 174, + 157, + 175 + ], + [ + 137, + 156, + 155 + ], + [ + 156, + 173, + 155 + ], + [ + 156, + 137, + 138 + ], + [ + 156, + 174, + 173 + ], + [ + 157, + 156, + 138 + ], + [ + 174, + 156, + 157 + ], + [ + 79, + 78, + 60 + ], + [ + 117, + 116, + 98 + ], + [ + 116, + 117, + 134 + ], + [ + 133, + 116, + 134 + ], + [ + 116, + 133, + 115 + ], + [ + 117, + 135, + 134 + ], + [ + 135, + 153, + 134 + ], + [ + 153, + 135, + 154 + ], + [ + 135, + 136, + 154 + ], + [ + 135, + 117, + 118 + ], + [ + 136, + 135, + 118 + ], + [ + 153, + 152, + 134 + ], + [ + 133, + 152, + 151 + ], + [ + 152, + 133, + 134 + ], + [ + 151, + 152, + 169 + ], + [ + 170, + 171, + 189 + ], + [ + 170, + 153, + 171 + ], + [ + 188, + 170, + 189 + ], + [ + 170, + 152, + 153 + ], + [ + 170, + 188, + 169 + ], + [ + 152, + 170, + 169 + ], + [ + 132, + 151, + 150 + ], + [ + 132, + 133, + 151 + ], + [ + 133, + 132, + 115 + ], + [ + 132, + 114, + 115 + ], + [ + 132, + 150, + 131 + ], + [ + 114, + 132, + 131 + ], + [ + 77, + 78, + 95 + ], + [ + 94, + 77, + 95 + ], + [ + 59, + 77, + 58 + ], + [ + 77, + 59, + 78 + ], + [ + 77, + 76, + 58 + ], + [ + 77, + 94, + 76 + ], + [ + 30, + 48, + 29 + ], + [ + 31, + 30, + 13 + ], + [ + 30, + 31, + 49 + ], + [ + 48, + 30, + 49 + ], + [ + 82, + 81, + 63 + ], + [ + 81, + 82, + 99 + ], + [ + 81, + 62, + 63 + ], + [ + 62, + 81, + 80 + ], + [ + 81, + 98, + 80 + ], + [ + 81, + 99, + 98 + ], + [ + 44, + 61, + 43 + ], + [ + 62, + 61, + 44 + ], + [ + 61, + 60, + 43 + ], + [ + 61, + 79, + 60 + ], + [ + 61, + 62, + 80 + ], + [ + 79, + 61, + 80 + ], + [ + 10, + 9, + 8 + ], + [ + 9, + 26, + 8 + ], + [ + 45, + 62, + 44 + ], + [ + 26, + 45, + 44 + ], + [ + 62, + 45, + 63 + ], + [ + 45, + 46, + 63 + ], + [ + 128, + 111, + 129 + ], + [ + 111, + 112, + 129 + ], + [ + 92, + 111, + 110 + ], + [ + 111, + 128, + 110 + ], + [ + 111, + 92, + 93 + ], + [ + 112, + 111, + 93 + ], + [ + 183, + 184, + 201 + ], + [ + 184, + 166, + 185 + ], + [ + 166, + 184, + 165 + ], + [ + 184, + 183, + 165 + ], + [ + 221, + 203, + 222 + ], + [ + 239, + 221, + 222 + ], + [ + 221, + 239, + 238 + ], + [ + 220, + 221, + 238 + ], + [ + 261, + 279, + 278 + ], + [ + 279, + 297, + 278 + ], + [ + 279, + 262, + 280 + ], + [ + 279, + 261, + 262 + ], + [ + 298, + 279, + 280 + ], + [ + 297, + 279, + 298 + ], + [ + 311, + 330, + 329 + ], + [ + 330, + 311, + 312 + ], + [ + 311, + 294, + 312 + ], + [ + 311, + 329, + 310 + ], + [ + 293, + 311, + 310 + ], + [ + 294, + 311, + 293 + ], + [ + 68, + 69, + 87 + ], + [ + 68, + 51, + 69 + ], + [ + 86, + 68, + 87 + ], + [ + 68, + 86, + 67 + ], + [ + 67, + 66, + 49 + ], + [ + 66, + 48, + 49 + ], + [ + 83, + 82, + 64 + ], + [ + 82, + 83, + 100 + ], + [ + 140, + 141, + 159 + ], + [ + 140, + 159, + 158 + ], + [ + 139, + 140, + 158 + ], + [ + 140, + 139, + 122 + ], + [ + 88, + 70, + 89 + ], + [ + 88, + 69, + 70 + ], + [ + 69, + 88, + 87 + ], + [ + 88, + 105, + 87 + ], + [ + 97, + 79, + 80 + ], + [ + 97, + 116, + 115 + ], + [ + 98, + 97, + 80 + ], + [ + 116, + 97, + 98 + ], + [ + 30, + 12, + 13 + ], + [ + 12, + 11, + 10 + ], + [ + 11, + 12, + 29 + ], + [ + 12, + 30, + 29 + ], + [ + 13, + 12, + 8 + ], + [ + 12, + 10, + 8 + ], + [ + 27, + 28, + 46 + ], + [ + 9, + 27, + 26 + ], + [ + 28, + 27, + 10 + ], + [ + 27, + 9, + 10 + ], + [ + 45, + 27, + 46 + ], + [ + 27, + 45, + 26 + ], + [ + 202, + 220, + 201 + ], + [ + 184, + 202, + 201 + ], + [ + 221, + 202, + 203 + ], + [ + 202, + 221, + 220 + ], + [ + 203, + 202, + 185 + ], + [ + 202, + 184, + 185 + ], + [ + 50, + 68, + 67 + ], + [ + 68, + 50, + 51 + ], + [ + 31, + 50, + 49 + ], + [ + 50, + 67, + 49 + ], + [ + 50, + 31, + 32 + ], + [ + 51, + 50, + 32 + ], + [ + 48, + 65, + 47 + ], + [ + 66, + 65, + 48 + ], + [ + 65, + 64, + 47 + ], + [ + 65, + 83, + 64 + ], + [ + 65, + 66, + 84 + ], + [ + 83, + 65, + 84 + ], + [ + 84, + 85, + 102 + ], + [ + 66, + 85, + 84 + ], + [ + 102, + 85, + 103 + ], + [ + 85, + 86, + 103 + ], + [ + 86, + 85, + 67 + ], + [ + 85, + 66, + 67 + ], + [ + 83, + 101, + 100 + ], + [ + 101, + 119, + 100 + ], + [ + 101, + 84, + 102 + ], + [ + 101, + 83, + 84 + ], + [ + 119, + 101, + 120 + ], + [ + 101, + 102, + 120 + ], + [ + 104, + 86, + 87 + ], + [ + 105, + 104, + 87 + ], + [ + 104, + 122, + 103 + ], + [ + 86, + 104, + 103 + ], + [ + 123, + 124, + 141 + ], + [ + 123, + 105, + 124 + ], + [ + 123, + 140, + 122 + ], + [ + 140, + 123, + 141 + ], + [ + 104, + 123, + 122 + ], + [ + 123, + 104, + 105 + ], + [ + 105, + 106, + 124 + ], + [ + 88, + 106, + 105 + ], + [ + 125, + 106, + 107 + ], + [ + 124, + 106, + 125 + ], + [ + 107, + 106, + 89 + ], + [ + 106, + 88, + 89 + ], + [ + 114, + 96, + 115 + ], + [ + 96, + 97, + 115 + ], + [ + 96, + 114, + 95 + ], + [ + 97, + 96, + 79 + ], + [ + 78, + 96, + 95 + ], + [ + 79, + 96, + 78 + ] + ], + "domains": { + "A": [ + 19, + 53, + 83, + 84, + 85, + 86, + 87, + 117, + 118, + 119, + 120, + 121, + 160, + 161, + 164, + 166, + 175, + 176, + 235, + 237, + 238, + 239, + 240, + 265, + 269, + 270, + 271, + 272, + 273, + 274, + 275, + 276, + 277, + 278, + 279, + 281, + 282, + 288, + 289, + 290, + 291, + 292, + 293, + 298, + 299, + 302, + 429, + 430, + 431, + 432, + 433, + 434, + 435, + 437, + 440, + 441, + 442, + 443, + 444, + 445, + 446, + 447, + 448, + 455, + 456, + 457, + 458, + 459, + 494, + 495, + 498, + 499, + 500, + 501, + 502, + 503, + 504, + 505, + 506, + 507, + 508, + 509, + 510, + 512, + 513, + 517, + 528, + 529, + 530, + 531, + 532, + 533, + 589, + 590, + 591, + 592, + 595, + 596, + 597, + 598, + 599, + 600, + 601, + 602, + 603, + 604, + 605, + 606, + 607, + 608, + 610, + 611, + 612, + 613, + 614, + 615, + 616, + 617, + 618, + 619, + 620, + 621, + 622, + 623, + 624, + 625, + 626, + 627, + 628, + 629, + 630, + 631, + 632, + 633, + 634, + 635, + 674, + 675, + 676, + 678, + 679, + 707, + 724, + 726, + 750, + 768 + ], + "B": [ + 6, + 7, + 36, + 54, + 55, + 56, + 57, + 58, + 59, + 60, + 61, + 62, + 63, + 64, + 65, + 66, + 67, + 68, + 69, + 70, + 71, + 72, + 88, + 89, + 90, + 124, + 127, + 128, + 159, + 162, + 163, + 165, + 167, + 177, + 179, + 180, + 181, + 182, + 183, + 184, + 185, + 186, + 187, + 188, + 189, + 190, + 191, + 192, + 193, + 229, + 230, + 231, + 232, + 233, + 234, + 280, + 283, + 284, + 285, + 286, + 287, + 294, + 295, + 296, + 297, + 300, + 301, + 303, + 310, + 311, + 312, + 313, + 314, + 315, + 316, + 317, + 318, + 319, + 320, + 321, + 322, + 335, + 336, + 337, + 338, + 339, + 340, + 345, + 350, + 351, + 352, + 353, + 354, + 355, + 356, + 357, + 358, + 359, + 360, + 361, + 362, + 363, + 364, + 365, + 366, + 367, + 368, + 369, + 370, + 371, + 372, + 373, + 374, + 375, + 376, + 377, + 378, + 379, + 380, + 381, + 382, + 383, + 384, + 385, + 386, + 387, + 388, + 389, + 390, + 391, + 392, + 393, + 394, + 395, + 396, + 397, + 398, + 399, + 400, + 401, + 402, + 403, + 404, + 405, + 406, + 407, + 408, + 425, + 426, + 427, + 428, + 511, + 514, + 515, + 516, + 518, + 519, + 520, + 521, + 522, + 523, + 524, + 525, + 526, + 527, + 534, + 535, + 536, + 537, + 538, + 539, + 540, + 541, + 542, + 543, + 544, + 545, + 546, + 547, + 548, + 549, + 550, + 551, + 552, + 553, + 554, + 555, + 556, + 557, + 558, + 559, + 560, + 561, + 562, + 563, + 564, + 565, + 566, + 567, + 568, + 569, + 570, + 571, + 572, + 573, + 574, + 575, + 576, + 577, + 578, + 579, + 580, + 581, + 680, + 681, + 682, + 683, + 684, + 685, + 686, + 687, + 688, + 689 + ], + "C": [ + 1, + 2, + 3, + 5, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 20, + 21, + 22, + 23, + 24, + 25, + 26, + 27, + 28, + 29, + 30, + 31, + 32, + 33, + 43, + 44, + 45, + 46, + 47, + 48, + 49, + 50, + 51, + 73, + 74, + 75, + 77, + 78, + 79, + 80, + 81, + 82, + 94, + 95, + 96, + 99, + 101, + 102, + 103, + 104, + 105, + 106, + 107, + 108, + 109, + 110, + 111, + 112, + 113, + 114, + 115, + 116, + 122, + 123, + 125, + 126, + 132, + 133, + 134, + 135, + 136, + 137, + 138, + 139, + 140, + 141, + 142, + 143, + 144, + 145, + 146, + 147, + 148, + 149, + 150, + 151, + 152, + 153, + 154, + 155, + 156, + 157, + 158, + 194, + 195, + 196, + 197, + 198, + 199, + 200, + 201, + 202, + 203, + 204, + 205, + 206, + 217, + 218, + 219, + 220, + 221, + 222, + 223, + 224, + 225, + 226, + 227, + 228, + 246, + 247, + 248, + 249, + 250, + 251, + 252, + 253, + 254, + 255, + 256, + 257, + 258, + 259, + 260, + 261, + 262, + 263, + 264, + 266, + 267, + 268, + 304, + 305, + 306, + 307, + 308, + 309, + 323, + 324, + 325, + 326, + 327, + 328, + 329, + 330, + 331, + 332, + 333, + 334, + 409, + 410, + 411, + 412, + 413, + 414, + 419, + 420, + 421, + 422, + 423, + 424, + 461, + 475, + 476, + 477, + 478, + 479, + 480, + 481, + 482, + 483, + 484, + 485, + 486, + 487, + 488, + 489, + 490, + 491, + 492, + 493, + 496, + 497, + 640, + 641, + 664, + 665, + 666, + 667, + 668, + 669, + 670, + 671, + 672, + 673, + 677, + 722, + 723, + 725, + 727 + ], + "D": [ + 0, + 4, + 34, + 35, + 37, + 38, + 39, + 40, + 41, + 42, + 52, + 76, + 91, + 92, + 93, + 97, + 98, + 100, + 129, + 130, + 131, + 168, + 169, + 170, + 171, + 172, + 173, + 174, + 178, + 207, + 208, + 209, + 210, + 211, + 212, + 213, + 214, + 215, + 216, + 236, + 241, + 242, + 243, + 244, + 245, + 341, + 342, + 343, + 344, + 346, + 347, + 348, + 349, + 415, + 416, + 417, + 418, + 436, + 438, + 439, + 449, + 450, + 451, + 452, + 453, + 454, + 460, + 462, + 463, + 464, + 465, + 466, + 467, + 468, + 469, + 470, + 471, + 472, + 473, + 474, + 582, + 583, + 584, + 585, + 586, + 587, + 588, + 593, + 594, + 609, + 636, + 637, + 638, + 639, + 642, + 643, + 644, + 645, + 646, + 647, + 648, + 649, + 650, + 651, + 652, + 653, + 654, + 655, + 656, + 657, + 658, + 659, + 660, + 661, + 662, + 663, + 690, + 691, + 692, + 693, + 694, + 695, + 696, + 697, + 698, + 699, + 700, + 701, + 702, + 703, + 704, + 705, + 706, + 708, + 709, + 710, + 711, + 712, + 713, + 714, + 715, + 716, + 717, + 718, + 719, + 720, + 721, + 728, + 729, + 730, + 731, + 732, + 733, + 734, + 735, + 736, + 737, + 738, + 739, + 740, + 741, + 742, + 743, + 744, + 745, + 746, + 747, + 748, + 749, + 751, + 752, + 753, + 754, + 755, + 756, + 757, + 758, + 759, + 760, + 761, + 762, + 763, + 764, + 765, + 766, + 767, + 769, + 770, + 771, + 772, + 773 + ] + }, + "domain_centroids": { + "A": [ + -0.1029, + -0.1178 + ], + "B": [ + 0.6824, + 1.1166 + ], + "C": [ + -1.4938, + -0.0269 + ], + "D": [ + 0.5989, + -1.2335 + ] + }, + "centroid_a": [ + -0.1029, + -0.1178 + ], + "centroid_b": [ + 0.6824, + 1.1166 + ], + "test_point": [ + 0.3, + 0.3 + ], + "view": 1.8 +} \ No newline at end of file diff --git a/.claude/skills/cetz-figures/examples/domain-demo.png b/.claude/skills/cetz-figures/examples/domain-demo.png new file mode 100644 index 00000000..f439f354 Binary files /dev/null and b/.claude/skills/cetz-figures/examples/domain-demo.png differ diff --git a/.claude/skills/cetz-figures/examples/domain-demo.typ b/.claude/skills/cetz-figures/examples/domain-demo.typ new file mode 100644 index 00000000..52ef2c2b --- /dev/null +++ b/.claude/skills/cetz-figures/examples/domain-demo.typ @@ -0,0 +1,124 @@ +#import "@preview/cetz:0.3.4" + +#set page(width: auto, height: auto, margin: 8pt) +#set text(size: 9pt) + +// ── Mesh data ───────────────────────────────────────────────────────── +#let mesh = json("domain-demo-data.json") + +#let pt(i) = { + let v = mesh.vertices.at(i) + (v.at(0), v.at(1)) +} + +#let VIEW = mesh.view +#let ca = (mesh.centroid_a.at(0), mesh.centroid_a.at(1)) +#let cb = (mesh.centroid_b.at(0), mesh.centroid_b.at(1)) +#let cc = (mesh.domain_centroids.C.at(0), mesh.domain_centroids.C.at(1)) +#let cd = (mesh.domain_centroids.D.at(0), mesh.domain_centroids.D.at(1)) +#let xp = (mesh.test_point.at(0), mesh.test_point.at(1)) + +// ── Helpers ────────────────────────────────────────────────────────── +#let tri-centroid(tri) = { + let a = pt(tri.at(0)) + let b = pt(tri.at(1)) + let c = pt(tri.at(2)) + ((a.at(0) + b.at(0) + c.at(0)) / 3, + (a.at(1) + b.at(1) + c.at(1)) / 3) +} + +#let in-view(p) = { + calc.abs(p.at(0)) <= VIEW and calc.abs(p.at(1)) <= VIEW +} + +// ── Domain lookup ──────────────────────────────────────────────────── +#let domain-a-set = mesh.domains.A +#let domain-b-set = mesh.domains.B +#let domain-c-set = mesh.domains.C +#let domain-d-set = mesh.domains.D + +#let get-domain(idx) = { + if domain-a-set.contains(idx) { "A" } + else if domain-b-set.contains(idx) { "B" } + else if domain-c-set.contains(idx) { "C" } + else { "D" } +} + +// ── Colours ────────────────────────────────────────────────────────── +#let fill-a = rgb(70, 110, 180, 45) // blue tint +#let fill-b = rgb(180, 70, 70, 45) // red tint +#let fill-c = rgb(90, 155, 90, 40) // green tint +#let fill-d = rgb(180, 160, 70, 40) // yellow tint + +#let stroke-mesh = 0.3pt + rgb("#b0b0b0") + +#let colour-a = rgb("#2563eb") // blue +#let colour-b = rgb("#dc2626") // red +#let colour-c = rgb("#16a34a") // green +#let colour-d = rgb("#ca8a04") // amber +#let colour-xp = rgb("#7c3aed") // violet + +// ── Dot helper ─────────────────────────────────────────────────────── +#let dot(p, label: none, direction: (0.06, 0.06), align-to: "west", + radius: 0.028, colour: black, label-colour: none) = { + import cetz.draw: * + circle(p, radius: radius, fill: colour, stroke: none) + if label != none { + let lcolour = if label-colour == none { black } else { label-colour } + content( + (p.at(0) + direction.at(0), p.at(1) + direction.at(1)), + text(fill: lcolour, label), + anchor: align-to, + ) + } +} + +// ── Figure ─────────────────────────────────────────────────────────── +#box( + clip: true, + width: 10cm, + height: 10cm, + stroke: 0.5pt + luma(50%), + align(center + horizon, cetz.canvas(length: 2.5cm, { + import cetz.draw: * + + // 1. Draw all triangles with domain shading + for (idx, tri) in mesh.triangles.enumerate() { + let c = tri-centroid(tri) + if in-view(c) { + let d = get-domain(idx) + let fill-col = if d == "A" { fill-a } + else if d == "B" { fill-b } + else if d == "C" { fill-c } + else { fill-d } + line( + pt(tri.at(0)), pt(tri.at(1)), pt(tri.at(2)), + close: true, + fill: fill-col, + stroke: stroke-mesh, + ) + } + } + + // 2. Dashed lines from test point to A and B centroids + let dashed-a = (paint: colour-a, thickness: 0.8pt, dash: "dashed") + let dashed-b = (paint: colour-b, thickness: 0.8pt, dash: "dashed") + line(xp, ca, stroke: dashed-a) + line(xp, cb, stroke: dashed-b) + + // 3. All domain centroids + dot(ca, label: $c _ A$, direction: (-0.12, -0.08), align-to: "east", + radius: 0.032, colour: colour-a, label-colour: colour-a) + dot(cb, label: $c _ B$, direction: (0.08, 0.06), align-to: "west", + radius: 0.032, colour: colour-b, label-colour: colour-b) + dot(cc, label: $c _ C$, direction: (-0.10, 0.06), align-to: "east", + radius: 0.032, colour: colour-c, label-colour: colour-c) + dot(cd, label: $c _ D$, direction: (0.08, -0.08), align-to: "west", + radius: 0.032, colour: colour-d, label-colour: colour-d) + + // 4. Test particle + dot(xp, label: $x _ p$, direction: (0.08, 0.06), align-to: "west", + radius: 0.028, colour: colour-xp, label-colour: colour-xp) + + })), +) diff --git a/.claude/skills/cetz-figures/examples/facet-vs-true-normals.png b/.claude/skills/cetz-figures/examples/facet-vs-true-normals.png new file mode 100644 index 00000000..dd387c66 Binary files /dev/null and b/.claude/skills/cetz-figures/examples/facet-vs-true-normals.png differ diff --git a/.claude/skills/cetz-figures/examples/facet-vs-true-normals.typ b/.claude/skills/cetz-figures/examples/facet-vs-true-normals.typ new file mode 100644 index 00000000..8f1b7169 --- /dev/null +++ b/.claude/skills/cetz-figures/examples/facet-vs-true-normals.typ @@ -0,0 +1,119 @@ +#import "@preview/cetz:0.3.4" + +// Explicit page size + aligned-centred canvas lets us control padding +// independently of whatever bbox cetz reports for its content. +#set page(width: 12cm, height: 8cm, margin: 6pt) +#set text(size: 10pt) + +#let data = json("curved-bc-data.json") + +// ── Colours ─────────────────────────────────────────────────────────── +#let arc-colour = rgb("#b0b0b0") // true smooth curve +#let facet-colour = rgb("#1f3a6b") // navy — mesh facets +#let gamma-colour = rgb("#c2410c") // rust — PETSc facet normal +#let true-colour = rgb("#059669") // emerald — true surface normal +#let centre-colour = black + +// ── Geometry access ─────────────────────────────────────────────────── +#let centre = (data.centre.at(0), data.centre.at(1)) +#let radius-end = (data.radius_end.at(0), data.radius_end.at(1)) + +// ── Figure ──────────────────────────────────────────────────────────── +#align(center + horizon, cetz.canvas(length: 3cm, { + import cetz.draw: * + + // 1. True smooth arc — thin dashed grey polyline. + let arc-stroke = (paint: arc-colour, thickness: 0.6pt, dash: "dashed") + for i in range(data.arc_points.len() - 1) { + let p0 = data.arc_points.at(i) + let p1 = data.arc_points.at(i + 1) + line((p0.at(0), p0.at(1)), (p1.at(0), p1.at(1)), stroke: arc-stroke) + } + + // 2. Radius-of-curvature indicator from centre to the arc midpoint. + line(centre, radius-end, + stroke: (paint: arc-colour, thickness: 0.5pt, dash: "dotted")) + // Label the radius midway along it, offset perpendicular. + let r-label-pos = ( + 0.5 * (centre.at(0) + radius-end.at(0)) + 0.07, + 0.5 * (centre.at(1) + radius-end.at(1)), + ) + content(r-label-pos, text(fill: arc-colour, size: 10pt, $R$), + anchor: "west") + + // 3. Facet polyline — thin black. The control points (vertex dots, + // Gauss-point dots) carry the visual weight; the segments between + // them are just the domain boundary. + for i in range(data.facet_vertices.len() - 1) { + let p0 = data.facet_vertices.at(i) + let p1 = data.facet_vertices.at(i + 1) + line((p0.at(0), p0.at(1)), (p1.at(0), p1.at(1)), + stroke: 0.6pt + black) + } + // Facet vertex markers + for v in data.facet_vertices { + circle((v.at(0), v.at(1)), radius: 0.038, + fill: facet-colour, stroke: none) + } + + // 4. Quadrature points with their two normals. + // - true normal (emerald, dashed) ← what free-slip wants + // - facet normal (rust, solid) ← what mesh.Gamma gives + let arrow-len = 0.28 + let true-stroke = (paint: true-colour, thickness: 0.9pt, dash: "dashed") + let gamma-stroke = 1.3pt + gamma-colour + for q in data.quadrature { + let pos = (q.pos.at(0), q.pos.at(1)) + let tn = q.true_normal + let fn = q.facet_normal + + // Facet normal arrow drawn FIRST, so the dashed true normal + // layers on top. Where the two coincide (facet midpoints), the + // green dashes interleave with the rust underneath and both + // colours stay visible — it reads as "they overlap here". + let tip-f = (pos.at(0) + arrow-len * fn.at(0), + pos.at(1) + arrow-len * fn.at(1)) + line(pos, tip-f, stroke: gamma-stroke, + mark: (end: ">", fill: gamma-colour)) + + // True normal arrow on top + let tip-t = (pos.at(0) + arrow-len * tn.at(0), + pos.at(1) + arrow-len * tn.at(1)) + line(pos, tip-t, stroke: true-stroke, + mark: (end: ">", fill: true-colour)) + + // Small quadrature dot + circle(pos, radius: 0.024, fill: black, stroke: none) + } + + // 5. Centre dot + label + circle(centre, radius: 0.038, fill: centre-colour, stroke: none) + content((centre.at(0) + 0.06, centre.at(1) - 0.06), + $O$, anchor: "north-west") + + // 6. Legend in the empty space below the arc. Short sample strokes + // next to each label, so the figure is self-explaining without + // leader lines or per-arrow annotations. + let lx = 0.25 // legend x start (inside the open area) + let ly = 0.55 // top row y + let row = 0.22 // row spacing + let sample-len = 0.28 + + // Row 1: facet normal sample + line((lx, ly), (lx + sample-len, ly), + stroke: gamma-stroke, + mark: (end: ">", fill: gamma-colour)) + content((lx + sample-len + 0.08, ly), + text(fill: gamma-colour, size: 9.5pt, + $hat(n)_Gamma$ + [ (facet, `mesh.Gamma`)]), + anchor: "west") + + // Row 2: true normal sample + line((lx, ly - row), (lx + sample-len, ly - row), + stroke: true-stroke, + mark: (end: ">", fill: true-colour)) + content((lx + sample-len + 0.08, ly - row), + text(fill: true-colour, size: 9.5pt, + $hat(n)_"true"$ + [ (smooth surface)]), + anchor: "west") +})) diff --git a/.claude/skills/cetz-figures/examples/generate-curved-bc-data.py b/.claude/skills/cetz-figures/examples/generate-curved-bc-data.py new file mode 100644 index 00000000..5af372f8 --- /dev/null +++ b/.claude/skills/cetz-figures/examples/generate-curved-bc-data.py @@ -0,0 +1,127 @@ +""" +Data for the "facet normals vs. true normals" figure in +`docs/advanced/curved-boundary-conditions.md`. + +The figure illustrates why `mesh.Gamma` (PETSc facet normals) diverges +from the true smooth-surface normal on a curved boundary: + + - Three straight facets approximate a circular arc. + - At each Gauss quadrature point (2 per facet), the facet normal + is constant across the facet while the true (radial) normal + rotates with position. + +Output schema: + + { + "centre": [cx, cy], + "radius": R, + "arc_points": [[x, y], ...], # densely sampled true arc + "facet_vertices": [[x, y], ...], # points on the circle + "quadrature": [ + {"pos": [x, y], + "facet_normal": [nx, ny], + "true_normal": [nx, ny], + "facet_idx": int}, + ... + ], + "radius_end": [x, y] # arc point at mid-angle + # (end of the radius indicator) + } +""" +import json +import math +from pathlib import Path + +CENTRE = (0.0, 0.0) +RADIUS = 1.5 +ANGLE_START = 30.0 # degrees — chosen so the error angle at +ANGLE_END = 150.0 # quadrature points is visually obvious (~12°) +N_FACETS = 3 +ARC_SAMPLES = 120 + + +def point_at_angle(deg, r=RADIUS): + rad = math.radians(deg) + return (CENTRE[0] + r * math.cos(rad), + CENTRE[1] + r * math.sin(rad)) + + +# 1. Dense sampling of the true arc (for a dashed polyline in Typst) +arc_points = [ + list(point_at_angle( + ANGLE_START + (i / ARC_SAMPLES) * (ANGLE_END - ANGLE_START) + )) + for i in range(ARC_SAMPLES + 1) +] + +# 2. Facet vertices — N_FACETS + 1 points evenly spaced on the arc +facet_vertices = [ + list(point_at_angle( + ANGLE_START + (i / N_FACETS) * (ANGLE_END - ANGLE_START) + )) + for i in range(N_FACETS + 1) +] + +# 3. Three-point Gauss–Legendre quadrature on [-1, 1]: 0, ±sqrt(3/5). +# Mapped to t ∈ [0, 1]: 0.5 ± sqrt(3/5)/2 and 0.5. +# The middle node lies exactly at the chord midpoint, where the +# facet normal and the radial true normal coincide — that's the +# "error vanishes at facet midpoint" case the figure needs to show. +_HALF = 0.5 * math.sqrt(0.6) +GAUSS_T = ( + 0.5 - _HALF, + 0.5, + 0.5 + _HALF, +) + +quadrature = [] +for facet_idx in range(N_FACETS): + p0 = facet_vertices[facet_idx] + p1 = facet_vertices[facet_idx + 1] + dx, dy = p1[0] - p0[0], p1[1] - p0[1] + length = math.hypot(dx, dy) + + # Outward perpendicular to the chord: pick the rotation of (dx, dy) + # whose dot-product with (midpoint - centre) is positive. + mid = ((p0[0] + p1[0]) / 2, (p0[1] + p1[1]) / 2) + cand = (-dy / length, dx / length) + if (cand[0] * (mid[0] - CENTRE[0]) + + cand[1] * (mid[1] - CENTRE[1])) < 0: + cand = (dy / length, -dx / length) + facet_normal = cand + + for t in GAUSS_T: + pos = (p0[0] + t * dx, p0[1] + t * dy) + # True normal: unit radial vector from circle centre through pos. + rx, ry = pos[0] - CENTRE[0], pos[1] - CENTRE[1] + rlen = math.hypot(rx, ry) + true_normal = (rx / rlen, ry / rlen) + quadrature.append({ + "pos": [pos[0], pos[1]], + "facet_normal": [facet_normal[0], facet_normal[1]], + "true_normal": [true_normal[0], true_normal[1]], + "facet_idx": facet_idx, + }) + +# 4. Radius indicator: from centre to arc midpoint +radius_end = list(point_at_angle((ANGLE_START + ANGLE_END) / 2)) + +data = { + "centre": list(CENTRE), + "radius": RADIUS, + "arc_points": [[round(x, 4), round(y, 4)] for x, y in arc_points], + "facet_vertices": [[round(x, 4), round(y, 4)] for x, y in facet_vertices], + "quadrature": [ + {k: ([round(v[0], 4), round(v[1], 4)] if isinstance(v, list) else v) + for k, v in q.items()} + for q in quadrature + ], + "radius_end": [round(radius_end[0], 4), round(radius_end[1], 4)], +} + +out = Path(__file__).with_name("curved-bc-data.json") +out.write_text(json.dumps(data, indent=2)) +print(f"wrote {out.name}: " + f"{len(arc_points)} arc samples, " + f"{len(facet_vertices)} facet vertices, " + f"{len(quadrature)} quadrature points") diff --git a/.claude/skills/cetz-figures/examples/generate-domain-demo-data.py b/.claude/skills/cetz-figures/examples/generate-domain-demo-data.py new file mode 100644 index 00000000..1e4b9990 --- /dev/null +++ b/.claude/skills/cetz-figures/examples/generate-domain-demo-data.py @@ -0,0 +1,154 @@ +""" +Generate a 2D triangular mesh partitioned into two domains for the +parallel particle migration blog post figure. + +The mesh is larger than the element-level demo, with more equant +triangles. Two domains (A and B) are defined by a partition line. +A test particle is placed so it is closer to domain A's centroid +but actually inside domain B — illustrating why centroid-nearest +is insufficient for parallel migration. + +Output schema: + { + "vertices": [[x, y], ...], + "triangles": [[i, j, k], ...], + "domain_a": [tri_index, ...], # triangle indices in domain A + "domain_b": [tri_index, ...], # triangle indices in domain B + "centroid_a": [x, y], + "centroid_b": [x, y], + "test_point": [x, y] + } + +Run: python3 generate-domain-demo-data.py +""" +import json +from pathlib import Path + +import numpy as np +from scipy.spatial import Delaunay + +SEED = 42 +BOX = 2.5 +SPACING = 0.28 +JITTER = 0.10 + +rng = np.random.default_rng(SEED) + +# Jittered equilateral-triangular lattice +dy = SPACING * np.sqrt(3) / 2 +lattice = [] +ys = np.arange(-BOX, BOX + dy, dy) +for row, y in enumerate(ys): + offset = SPACING / 2 if row % 2 else 0.0 + for x in np.arange(-BOX + offset, BOX + 1e-9, SPACING): + lattice.append([x, y]) +lattice = np.array(lattice) + +# Jitter +lattice += rng.uniform(-JITTER * SPACING, JITTER * SPACING, size=lattice.shape) + +# Triangulate +tri = Delaunay(lattice) + +# Compute triangle centroids +centroids = lattice[tri.simplices].mean(axis=1) + +# Partition into four real domains, all with centroids: +# - Domain A (blue): L-shaped, left and lower +# - Domain B (red): upper-right plus top strip (contiguous) +# - Domain C (green): left strip +# - Domain D (yellow): bottom strip plus right edge of A's boundary + +def assign_domain(cx, cy): + """Assign a triangle centroid to a domain.""" + # Domain C (green): far left + if cx <= -1.2: + return "C" + # Domain D (yellow): far bottom, plus right boundary strip + if cy <= -1.2 or (cx > 1.4 and cy <= 0.1): + return "D" + # Domain B (red): upper-right quadrant plus entire top strip + if (cx > 0.1 and cy > 0.1) or cy > 1.4: + return "B" + # Domain A (blue): everything else (L-shaped, central-left and lower) + return "A" + +domains = {} +for label in ["A", "B", "C", "D"]: + domains[label] = [] + +for i, (cx, cy) in enumerate(centroids): + label = assign_domain(cx, cy) + domains[label].append(i) + +# Compute domain centroids +VIEW = 1.8 + +# Compute domain centroids using only triangles visible in the view, +# so the dots match what the reader sees in the clipped figure. +domain_centroids = {} +for label, indices in domains.items(): + visible = [i for i in indices + if abs(centroids[i, 0]) <= VIEW and abs(centroids[i, 1]) <= VIEW] + if visible: + domain_centroids[label] = centroids[visible].mean(axis=0) + elif indices: + domain_centroids[label] = centroids[indices].mean(axis=0) + +centroid_a = domain_centroids["A"] +centroid_b = domain_centroids["B"] + +# Place test particle: inside domain B but closer to centroid A. +test_point = [0.30, 0.30] + +dist_a = np.hypot(test_point[0] - centroid_a[0], test_point[1] - centroid_a[1]) +dist_b = np.hypot(test_point[0] - centroid_b[0], test_point[1] - centroid_b[1]) + +print(f"Domain A centroid: ({centroid_a[0]:.3f}, {centroid_a[1]:.3f})") +print(f"Domain B centroid: ({centroid_b[0]:.3f}, {centroid_b[1]:.3f})") +print(f"Test point: ({test_point[0]:.3f}, {test_point[1]:.3f})") +print(f"Distance to A: {dist_a:.3f}, Distance to B: {dist_b:.3f}") +print(f"Closer to A: {dist_a < dist_b}") + +if dist_a >= dist_b: + for ty in np.arange(0.12, 0.8, 0.02): + for tx in np.arange(0.12, 0.8, 0.02): + if assign_domain(tx, ty) != "B": + continue + da = np.hypot(tx - centroid_a[0], ty - centroid_a[1]) + db = np.hypot(tx - centroid_b[0], ty - centroid_b[1]) + if da < db: + test_point = [tx, ty] + print(f"Adjusted to ({tx:.3f}, {ty:.3f}), dist_a={da:.3f}, dist_b={db:.3f}") + break + else: + continue + break + +# All domain centroids +domain_centroids_out = {} +for label, indices in domains.items(): + if indices: + c = domain_centroids[label] + domain_centroids_out[label] = [round(float(c[0]), 4), round(float(c[1]), 4)] + +data = { + "vertices": [[round(float(x), 4), round(float(y), 4)] for x, y in lattice], + "triangles": [[int(i), int(j), int(k)] for i, j, k in tri.simplices], + "domains": {label: [int(i) for i in indices] for label, indices in domains.items()}, + "domain_centroids": domain_centroids_out, + "centroid_a": [round(float(centroid_a[0]), 4), round(float(centroid_a[1]), 4)], + "centroid_b": [round(float(centroid_b[0]), 4), round(float(centroid_b[1]), 4)], + "test_point": [round(float(test_point[0]), 4), round(float(test_point[1]), 4)], + "view": VIEW, +} + +out = Path(__file__).with_name("domain-demo-data.json") +out.write_text(json.dumps(data, indent=2)) + +print(f"\nWrote {out.name}: " + f"{len(data['vertices'])} vertices, " + f"{len(data['triangles'])} triangles") +for label, indices in domains.items(): + c = domain_centroids[label] + print(f" Domain {label}: {len(indices)} triangles, centroid ({c[0]:.3f}, {c[1]:.3f})") diff --git a/.claude/skills/cetz-figures/examples/generate-mesh-data.py b/.claude/skills/cetz-figures/examples/generate-mesh-data.py new file mode 100644 index 00000000..9fe4804a --- /dev/null +++ b/.claude/skills/cetz-figures/examples/generate-mesh-data.py @@ -0,0 +1,107 @@ +""" +Generate a deterministic 2D triangular mesh for the cetz demo figure. + +The central "highlight" triangle is forced into the Delaunay output by +seeding its vertices first and rejecting any other interior point that +would fall inside its circumcircle (Delaunay-empty-circle property). + +Output schema (identical to what a future underworld3 export will emit): + + { + "vertices": [[x, y], ...], # N x 2 + "triangles": [[i, j, k], ...], # M x 3 indices + "highlight": int # index into triangles + } + +Run: python3 generate-mesh-data.py +""" +import json +from pathlib import Path + +import numpy as np +from scipy.spatial import Delaunay + +SEED = 3 +BOX = 2.0 # point cloud extent — further out to ensure + # enough lattice points survive for cut edges on + # every side at the coarse spacing below. +SPACING = 1.25 # target edge length of the FE-like lattice — + # matched to the central triangle's average edge + # (~1.29) so it reads as a typical cell. +JITTER = 0.15 # fraction of SPACING — small randomness so the + # mesh reads as a real FE discretisation rather + # than a perfect lattice + +# Central triangle — deliberately asymmetric. +TRI = np.array([ + [-0.60, -0.38], + [ 0.58, -0.45], # v2 pulled in slightly from (0.62, -0.48) + # so it isn't flush with the frame. + [ 0.32, 0.70], # v3 shifted up-right so CA is the longest edge. +]) + + +def circumcircle(a, b, c): + """Return (cx, cy, r) of the circumscribed circle of triangle abc.""" + ax, ay = a; bx, by = b; cx_, cy_ = c + d = 2 * (ax * (by - cy_) + bx * (cy_ - ay) + cx_ * (ay - by)) + ux = ((ax**2 + ay**2) * (by - cy_) + + (bx**2 + by**2) * (cy_ - ay) + + (cx_**2 + cy_**2) * (ay - by)) / d + uy = ((ax**2 + ay**2) * (cx_ - bx) + + (bx**2 + by**2) * (ax - cx_) + + (cx_**2 + cy_**2) * (bx - ax)) / d + r = float(np.hypot(ux - ax, uy - ay)) + return float(ux), float(uy), r + + +cx, cy, r = circumcircle(*TRI) +SAFETY = 0.05 # keep candidates a little further out than the circle + +rng = np.random.default_rng(SEED) + +# Jittered equilateral-triangular lattice — visually uniform triangle +# sizes, like an FE mesh from gmsh with a single size field. Every other +# row is offset by half a spacing to build equilateral triangles. +dy = SPACING * np.sqrt(3) / 2 +lattice = [] +ys = np.arange(-BOX, BOX + dy, dy) +for row, y in enumerate(ys): + offset = SPACING / 2 if row % 2 else 0.0 + for x in np.arange(-BOX + offset, BOX + 1e-9, SPACING): + lattice.append([x, y]) +lattice = np.array(lattice) + +# Add a small jitter (after placing) so the mesh doesn't read as a +# perfect lattice — closer to a real FE mesh in appearance. +lattice += rng.uniform(-JITTER * SPACING, JITTER * SPACING, + size=lattice.shape) + +# Keep only lattice points outside the central triangle's circumcircle. +dist = np.hypot(lattice[:, 0] - cx, lattice[:, 1] - cy) +lattice = lattice[dist > r + SAFETY] + +points = np.vstack([TRI, lattice]) +tri = Delaunay(points) + +# Locate the central triangle among the simplices (any vertex order), +# then canonicalise its row to (0, 1, 2) so the Typst figure can rely +# on triangle[highlight].at(k) referring to the kth original TRI vertex. +target = {0, 1, 2} +highlight = next(i for i, s in enumerate(tri.simplices) if set(s) == target) +simplices = tri.simplices.copy() +simplices[highlight] = [0, 1, 2] + +data = { + "vertices": [[round(float(x), 4), round(float(y), 4)] for x, y in points], + "triangles": [[int(i), int(j), int(k)] for i, j, k in simplices], + "highlight": int(highlight), +} + +out = Path(__file__).with_name("mesh-data.json") +out.write_text(json.dumps(data, indent=2)) + +print(f"wrote {out.name}: " + f"{len(data['vertices'])} vertices, " + f"{len(data['triangles'])} triangles, " + f"highlight={highlight}") diff --git a/.claude/skills/cetz-figures/examples/mesh-data.json b/.claude/skills/cetz-figures/examples/mesh-data.json new file mode 100644 index 00000000..83d31e85 --- /dev/null +++ b/.claude/skills/cetz-figures/examples/mesh-data.json @@ -0,0 +1,252 @@ +{ + "vertices": [ + [ + -0.6, + -0.38 + ], + [ + 0.58, + -0.45 + ], + [ + 0.32, + 0.7 + ], + [ + -2.1554, + -2.0987 + ], + [ + -0.637, + -1.9692 + ], + [ + 0.3478, + -2.0251 + ], + [ + 1.7421, + -2.1276 + ], + [ + -1.287, + -1.0623 + ], + [ + -0.1658, + -0.9112 + ], + [ + 1.099, + -0.8849 + ], + [ + -1.9108, + 0.3362 + ], + [ + -0.8309, + 0.2208 + ], + [ + 1.5631, + 0.3426 + ], + [ + -1.4506, + 1.1778 + ], + [ + 0.0219, + 1.2795 + ], + [ + 1.1142, + 1.3501 + ], + [ + -2.1761, + 2.4077 + ], + [ + -0.7972, + 2.1767 + ], + [ + 0.5602, + 2.4919 + ], + [ + 1.6402, + 2.3789 + ] + ], + "triangles": [ + [ + 15, + 19, + 18 + ], + [ + 3, + 10, + 16 + ], + [ + 10, + 13, + 16 + ], + [ + 4, + 7, + 3 + ], + [ + 7, + 10, + 3 + ], + [ + 7, + 8, + 0 + ], + [ + 8, + 7, + 4 + ], + [ + 17, + 18, + 16 + ], + [ + 13, + 17, + 16 + ], + [ + 14, + 17, + 13 + ], + [ + 14, + 15, + 18 + ], + [ + 17, + 14, + 18 + ], + [ + 19, + 12, + 6 + ], + [ + 15, + 12, + 19 + ], + [ + 5, + 8, + 4 + ], + [ + 5, + 3, + 6 + ], + [ + 5, + 4, + 3 + ], + [ + 11, + 14, + 13 + ], + [ + 10, + 11, + 13 + ], + [ + 11, + 7, + 0 + ], + [ + 7, + 11, + 10 + ], + [ + 12, + 9, + 6 + ], + [ + 9, + 5, + 6 + ], + [ + 5, + 9, + 8 + ], + [ + 2, + 11, + 0 + ], + [ + 11, + 2, + 14 + ], + [ + 14, + 2, + 15 + ], + [ + 2, + 12, + 15 + ], + [ + 8, + 1, + 0 + ], + [ + 0, + 1, + 2 + ], + [ + 9, + 1, + 8 + ], + [ + 1, + 9, + 12 + ], + [ + 2, + 1, + 12 + ] + ], + "highlight": 29 +} \ No newline at end of file diff --git a/.claude/skills/cetz-figures/examples/mesh-demo.png b/.claude/skills/cetz-figures/examples/mesh-demo.png new file mode 100644 index 00000000..d634ae21 Binary files /dev/null and b/.claude/skills/cetz-figures/examples/mesh-demo.png differ diff --git a/.claude/skills/cetz-figures/examples/mesh-demo.typ b/.claude/skills/cetz-figures/examples/mesh-demo.typ new file mode 100644 index 00000000..331bdea8 --- /dev/null +++ b/.claude/skills/cetz-figures/examples/mesh-demo.typ @@ -0,0 +1,204 @@ +#import "@preview/cetz:0.3.4" + +#set page(width: auto, height: auto, margin: 8pt) +#set text(size: 10pt) + +// ── Mesh data ───────────────────────────────────────────────────────── +#let mesh = json("mesh-data.json") + +#let pt(i) = { + let v = mesh.vertices.at(i) + (v.at(0), v.at(1)) +} + +// The central triangle is a genuine Delaunay cell. Python forces it +// into the triangulation and canonicalises its vertex order so +// triangle[highlight].at(k) is always TRI[k]. +#let h = mesh.triangles.at(mesh.highlight) +#let a = pt(h.at(0)) +#let b = pt(h.at(1)) +#let c = pt(h.at(2)) + +#let centroid = ( + (a.at(0) + b.at(0) + c.at(0)) / 3, + (a.at(1) + b.at(1) + c.at(1)) / 3, +) + +// ── Vector helpers ──────────────────────────────────────────────────── +#let vsub(p, q) = (p.at(0) - q.at(0), p.at(1) - q.at(1)) +#let vadd(p, q) = (p.at(0) + q.at(0), p.at(1) + q.at(1)) +#let vscale(p, s) = (s * p.at(0), s * p.at(1)) +#let vlen(p) = calc.sqrt(p.at(0) * p.at(0) + p.at(1) * p.at(1)) +#let vnorm(p) = { let m = vlen(p); (p.at(0) / m, p.at(1) / m) } +#let vmid(p, q) = ((p.at(0) + q.at(0)) / 2, (p.at(1) + q.at(1)) / 2) +#let dist(p, q) = vlen(vsub(p, q)) + +// ── Edge-midpoint markers (pair per edge: inside + outside) ─────────── +// Dot radius is 0.022 world units; gap 0.030 places the dot just clear. +#let edge-marker(p, q, gap: 0.030, side: "inside") = { + let m = vmid(p, q) + let dir = vnorm(vsub(centroid, m)) + let sign = if side == "outside" { -1.0 } else { 1.0 } + vadd(m, vscale(dir, sign * gap)) +} + +#let m-ab = edge-marker(a, b) +#let m-bc = edge-marker(b, c) +#let m-ca = edge-marker(c, a) + +#let m-ab-out = edge-marker(a, b, side: "outside") +#let m-bc-out = edge-marker(b, c, side: "outside") +#let m-ca-out = edge-marker(c, a, side: "outside") + +// For each face we store (endpoint0, endpoint1, inside-marker, outside-marker). +#let faces = ( + (a, b, m-ab, m-ab-out), + (b, c, m-bc, m-bc-out), + (c, a, m-ca, m-ca-out), +) + +// 2D cross product — used for the side-of-line test below. +#let cross2d(u, v) = u.at(0) * v.at(1) - u.at(1) * v.at(0) + +// Per-face rule: connect the test point to the *inside* marker if the +// test point lies on the same side of the edge as the centroid, else +// to the *outside* marker. A point inside the triangle lands on the +// centroid side of all three edges (→ all black); a point outside +// crosses at least one edge to the other side (→ at least one red). +#let marker-for(test, e0, e1, inside-marker, outside-marker) = { + let edge = vsub(e1, e0) + let s-test = cross2d(edge, vsub(test, e0)) + let s-ref = cross2d(edge, vsub(centroid, e0)) + if s-test * s-ref > 0 { inside-marker } else { outside-marker } +} + +// ── Sample (test) points ────────────────────────────────────────────── +// x_q: inside the triangle, pulled toward v2 so it sits closer to AB +// and BC, but shifted left enough that its line to the CA midpoint +// clears the centroid. +#let x-q = ( + 0.25 * a.at(0) + 0.60 * b.at(0) + 0.15 * c.at(0), + 0.25 * a.at(1) + 0.60 * b.at(1) + 0.15 * c.at(1), +) +// x_p: outside CA, nudged right so the connector to the AB midpoint +// doesn't graze the CA outside marker en route. +#let x-p = (-0.15, 0.80) + +// ── Line-through-edge extended to the viewport boundary ─────────────── +// Returns (before, after) — the two points where the infinite line +// through (p, q) crosses the axis-aligned box [-VIEW, VIEW]^2. +#let VIEW = 1.0 +#let extend-to-viewport(p, q) = { + let dx = q.at(0) - p.at(0) + let dy = q.at(1) - p.at(1) + let BIG = 1e6 + let tx-lo = if calc.abs(dx) < 1e-9 { -BIG } else { (-VIEW - p.at(0)) / dx } + let tx-hi = if calc.abs(dx) < 1e-9 { BIG } else { ( VIEW - p.at(0)) / dx } + let ty-lo = if calc.abs(dy) < 1e-9 { -BIG } else { (-VIEW - p.at(1)) / dy } + let ty-hi = if calc.abs(dy) < 1e-9 { BIG } else { ( VIEW - p.at(1)) / dy } + let tx-min = calc.min(tx-lo, tx-hi) + let tx-max = calc.max(tx-lo, tx-hi) + let ty-min = calc.min(ty-lo, ty-hi) + let ty-max = calc.max(ty-lo, ty-hi) + let t-enter = calc.max(tx-min, ty-min) + let t-exit = calc.min(tx-max, ty-max) + ( + (p.at(0) + t-enter * dx, p.at(1) + t-enter * dy), + (p.at(0) + t-exit * dx, p.at(1) + t-exit * dy), + ) +} + +// ── Colours ─────────────────────────────────────────────────────────── +#let tri-stroke = rgb("#1f3a6b") // navy — triangle +#let tri-fill = rgb(70, 120, 180, 55) // tinted navy +#let out-marker = rgb("#c2410c") // rust — outside midpoint dots +#let x-q-colour = rgb("#059669") // emerald — inside sample +#let x-p-colour = rgb("#7c3aed") // violet — outside sample + +// ── Labelled dot helper ─────────────────────────────────────────────── +#let dot(p, label: none, direction: (0.08, 0.08), align-to: "west", + radius: 0.022, colour: black, label-colour: none) = { + import cetz.draw: * + circle(p, radius: radius, fill: colour, stroke: none) + if label != none { + let lcolour = if label-colour == none { black } else { label-colour } + content( + (p.at(0) + direction.at(0), p.at(1) + direction.at(1)), + text(fill: lcolour, label), + anchor: align-to, + ) + } +} + +#let outward(p, gap: 0.11) = vscale(vnorm(vsub(p, centroid)), gap) + +// ── Figure ──────────────────────────────────────────────────────────── +#box( + clip: true, + width: 8cm, + height: 8cm, + stroke: 0.5pt + luma(50%), + align(center + horizon, cetz.canvas(length: 4cm, { + import cetz.draw: * + + // 1. Background mesh + for tri in mesh.triangles { + line( + pt(tri.at(0)), pt(tri.at(1)), pt(tri.at(2)), + close: true, + stroke: 0.4pt + rgb("#c0c0c0"), + ) + } + + // 2. Dotted extensions — triangle edges continued to the viewport. + let dotted = (paint: tri-stroke, thickness: 0.5pt, dash: "dotted") + for pair in ((a, b), (b, c), (c, a)) { + let p = pair.at(0) + let q = pair.at(1) + let ends = extend-to-viewport(p, q) + line(ends.at(0), p, stroke: dotted) + line(q, ends.at(1), stroke: dotted) + } + + // 3. Main triangle — fill + bold stroke on top of mesh and dots. + line(a, b, c, close: true, fill: tri-fill, stroke: 1.3pt + tri-stroke) + + // 4. Per-face connectors. Each test point draws three dashed lines, + // one per face, landing on black or red by the side-of-edge rule. + // All-black ⇒ inside; any-red ⇒ outside. + let dashed-stroke(col) = (paint: col, thickness: 0.7pt, dash: "dashed") + for face in faces { + let e0 = face.at(0) + let e1 = face.at(1) + let m-in = face.at(2) + let m-out = face.at(3) + line(x-q, marker-for(x-q, e0, e1, m-in, m-out), + stroke: dashed-stroke(x-q-colour)) + line(x-p, marker-for(x-p, e0, e1, m-in, m-out), + stroke: dashed-stroke(x-p-colour)) + } + + // 5. Vertex dots + labels outward from centroid. + let side-of(p) = if p.at(0) < centroid.at(0) { "east" } else { "west" } + dot(a, label: $v_1$, direction: outward(a), align-to: side-of(a)) + dot(b, label: $v_2$, direction: outward(b), align-to: side-of(b)) + dot(c, label: $v_3$, direction: outward(c), align-to: side-of(c)) + + // 6. Midpoint markers — inside (black) + outside companions (rust). + dot(m-ab); dot(m-bc); dot(m-ca) + dot(m-ab-out, colour: out-marker) + dot(m-bc-out, colour: out-marker) + dot(m-ca-out, colour: out-marker) + + // 7. Centroid + two sample points (each in its own colour). + // Centroid dot is smaller and its label tight against it so it + // reads as a derived reference rather than competing with the + // test points. + dot(centroid, label: $c$, direction: (0.055, 0), align-to: "west", + radius: 0.014) + dot(x-q, label: $x_q$, direction: (0.10, 0.05), align-to: "west", + colour: x-q-colour, label-colour: x-q-colour) + dot(x-p, label: $x_p$, direction: (0.10, 0.05), align-to: "west", + colour: x-p-colour, label-colour: x-p-colour) + })), +) diff --git a/.claude/skills/cetz-figures/underworld-bridge.md b/.claude/skills/cetz-figures/underworld-bridge.md new file mode 100644 index 00000000..f26c5791 --- /dev/null +++ b/.claude/skills/cetz-figures/underworld-bridge.md @@ -0,0 +1,75 @@ +# Bridging cetz figures to underworld3 mesh objects (sketch) + +This is a **sketch, not an implementation**. When we wire figures up to +real `underworld3` data, this is the contract the drawing code expects. + +## The JSON schema + +Any cetz figure that draws a mesh reads the same shape: + +```json +{ + "vertices": [[x, y], ...], // N x 2 floats + "triangles": [[i, j, k], ...], // M x 3 int indices into vertices + "highlight": 17 // optional: triangle of interest +} +``` + +Keep it minimal — cetz script should not have to do geometry or case +analysis. Anything derived (centroids, boundary detection, highlight +selection) is computed upstream and baked into the JSON. + +## How this maps to underworld3 + +Provisional — needs verification against the current `underworld3` mesh +API once we start using it: + +| JSON field | `underworld3` source | Notes | +|---|---|---| +| `vertices` | `mesh.data` (or `mesh.X.coords` — see repo style guide) | 2D only for this skill; 3D needs a projection | +| `triangles` | PETSc `DMPlex` cell-vertex connectivity | Simplex meshes only; quads/hexes need different rendering | +| `highlight` | Caller-provided or heuristic (e.g. closest-to-origin) | Chosen offline, not in Typst | + +Repo's data-access pattern guide: +`docs/developer/UW3_Style_and_Patterns_Guide.md`. + +## Python exporter sketch + +```python +# Not the current implementation — a target shape +def mesh_to_figure_json(mesh, *, highlight=None, path): + """Dump a uw.meshing object to the schema the cetz figures expect.""" + import json + # 2D coordinates + coords = mesh.X.coords # (N, 2) numpy, already preferred over mesh.data + cells = mesh.cell_vertex_indices # placeholder — actual accessor TBD + data = { + "vertices": [[round(float(x), 4), round(float(y), 4)] for x, y in coords], + "triangles": [[int(i) for i in cell] for cell in cells], + } + if highlight is not None: + data["highlight"] = int(highlight) + with open(path, "w") as f: + json.dump(data, f) +``` + +## What this intentionally does NOT do + +- **No coordinate transform / projection.** Caller passes whatever world + coordinates it wants; Typst draws them as-is. Transforming a spherical + mesh to a 2D projection is a Python concern. +- **No per-cell data.** This schema is for topology/geometry only. For + scalar fields, vector fields, colormaps → go the Python-generated-SVG + route and include via `#image()`. +- **No parallel-safe collection logic.** Caller is expected to gather the + mesh to rank 0 before exporting. Don't invent a parallel JSON writer. + +## Next steps (when we get here) + +1. Implement `mesh_to_figure_json` using the correct underworld3 accessors. + Verify which is idiomatic in the current codebase — `mesh.data` is + deprecated per `CLAUDE.md`. +2. Write a test figure consuming output from a small real mesh (e.g. one + of the existing example meshes). +3. Decide whether to ship the exporter as a `uw.utilities.export_figure` + helper or leave it as a per-paper script. diff --git a/.claude/skills/free-surface-convection/SKILL.md b/.claude/skills/free-surface-convection/SKILL.md new file mode 100644 index 00000000..e6687ef0 --- /dev/null +++ b/.claude/skills/free-surface-convection/SKILL.md @@ -0,0 +1,177 @@ +--- +name: free-surface-convection +description: The Underworld3 free-surface convection method we are hardening — the THREE-NUMBER pointwise topography integrator (held-lid stress equilibrium h_∞ + L-stable exponential relaxation), NOT FSSA. Reach for THIS before touching any free-surface / dynamic-topography convection run, choosing a surface-update scheme, or "stabilising" a free surface. It records the method, why FSSA is explicitly rejected, and the failure modes. +--- + +# free-surface-convection + +The free-surface scheme used in `~/+Simulations/FreeSurface/convection/fs4_compare.py` +and the design doc `docs/developer/design/FREESLIP_DYNAMIC_TOPOGRAPHY_FREESURFACE.md`. +**This is the method we are HARDENING — do not replace it; do not add FSSA.** + +> The 3-number integrator is necessary but NOT sufficient — see **Hardening strategies +> (2026-06)** below for the material-surface advection, tangential topography term, +> free-slip-inner nullspace, and graded/higher-order-mesh fixes that make it actually +> work. Reference impl + diagnostic tools live in `~/+Simulations/FreeSurface/convection/`. + +## ⚠️ NOT FSSA + +The `docs/examples/free_surface/advanced/Annulus*FS.py` examples use **FSSA** +(`add_natural_bc(δt·(Γ·v)Γ/2, "Upper")`). **That is NOT our method.** FSSA buys +stability by adding an implicit surface traction that **UNDER-deforms the surface** +— it trades accuracy for stability. Our scheme is designed to be stable **and** +accurate. If you find yourself adding `FSSA`, an `add_natural_bc` traction on the +free surface, or a `Gamma.dot(v)` stabiliser — STOP, you have the wrong method. +(Those example files are a template for a *different* approach, not this one.) + +## The three-number pointwise integrator (THE method) + +Two Stokes solves per step on the SAME mesh, then a pointwise surface update: + +1. **Free solve** — stress-free top (NO velocity BC on `Upper`; pressure datum is + pinned by the stress-free condition → no pressure nullspace). The surface + normal velocity `u_n` of this solve IS the kinematic rate `ḣ`. +2. **Held-lid solve** — a second Stokes solve with a RIGID free-slip held lid + (`u_n = 0`, via `add_nitsche_bc(Upper, local_h=True)` — see [[project_nitsche_local_h_pr275]]) + and a DRIVING-ONLY body force. Its surface normal stress `σ_nn` gives the + equilibrium topography `h_∞ = -(σ_nn - mean)/ρg`. (The free solve forces + `σ_nn = 0`, so the equilibrium MUST come from the held-lid stress.) +3. **Pointwise exp step**, per surface node, from THREE numbers (`h`, `ḣ=u_n`, + `h_∞`): + ``` + γ = ḣ / (h_∞ − h) # local relaxation rate, clamp γ ≥ 0 + h ← h_∞ + (h − h_∞)·exp(−γ·dt) + ``` + L-stable: the step is bounded between `h` and `h_∞`, so it **cannot overshoot** + regardless of a noisy local `γ` (no "drunken sailor"). 1 extra solve/step; + beats RK4 at large dt. NO per-node freeze-clamp (that was the old `relax` bug). +4. The nodal surface increment is **carried inward by a Laplacian diffuser** + (smooth, minimal mesh deformation — NOT full mmpde adaptation), then + `mesh.deform()`. Uniform meshes are fine; adaptivity is NOT required. + +Reference impl: `fs4_compare.py` → `_surface_step`, `_h_inf` (held-lid σ_nn via a +`Projection`), `_surf_un`, `_carry_diffuser`. The free-slip RIGID-top run (no +surface motion) is the reference; the free-surface run is the same driving solve +PLUS this surface update. See [[project_fs4_adaptive_2x2]], +[[project_stress_equilibrium_freesurface]], [[project_freeslip_topo_freesurface]]. + +## Performance + +- Free-slip (rigid top) stagnant-lid runs are FAST. The free-surface cost is the + extra held-lid solve **plus** that moving the surface forces a COLD-START Stokes + each step (can't warm-start across a deformed mesh). +- Use **uniform** meshes for this problem — the diffuser gives minimal deformation, + no mmpde needed. (FMG works on a uniform `refinement=N` hierarchy; scalar solvers + must avoid FMG — PETSc err62, issue underworldcode/underworld3#276.) + +## Hardening strategies (2026-06) — the integrator alone is not enough + +The 3-number integrator moves the surface correctly, but several *other* things must +be right or it runs away / tangles. All implemented in `fs4_compare.py` (flags noted). + +### 1. Material-surface advection — THE key fix (`--advect-velocity`) +The runaway (`u_n` 42→125→285→445, cold lid leaking in, plumes punching through) was +NOT an `h_∞`/BC bug (`h_∞` is verified correct, even in the stagnant FK lid — held-lid +free-slip is the EASY case there). The bug: the surface moves by the L-stable relaxed +rate `ũ_n = Δh/Δt ≤ u_n`, but T was advected with the **stress-free solve velocity** +(surface-normal = full `u_n`). Net material then crosses the surface. A free surface is +a MATERIAL boundary: advect T with a velocity whose surface-normal = `ũ_n`. Modes: +- `consistent` (the right way): a THIRD Stokes solve, same buoyancy, `v·n̂ = ũ_n` + PRESCRIBED at the surface (penalty), tangential stress-free. `ũ_n = (shape_new−shape0)/dt` + = the full ∂h/∂t at fixed θ (correct ALE target). +- `blend`: `α·v_free + (1−α)·v_held`, `α = φ1(γΔt) = (1−e^{−γΔt})/(γΔt)` (the exp-decay + time-average). By Stokes LINEARITY this *is* the prescribed-`ũ_n` solve for UNIFORM α + (and free for FK, which is linear in v). BUT the single mean-α collapse is NOT close + enough once γ varies per surface node — the planform diverges (mode-3 vs mode-1/2), + throughflow ~23 vs ~0.08. Per-node α breaks div-free (∇α·(v_free−v_held)). So the + per-node `consistent` 3rd solve is REQUIRED for structured planforms. +- `free`: advect with stress-free v (the inconsistent baseline — the runaway). + +### 2. Tangential topography advection (`--no-tangent-advect` to disable; default ON) +The pointwise relaxation omits the `v_t·∂_s h` term — a surface rotation/convergence +should carry the topography pattern along the surface; without it you get edge artefacts +where ∂_s h is large (plume-bulge edges). Fix = operator split per step: (1) departure- +point semi-Lagrangian transport of the surface shape in θ by `ω = v_t/r`, then (2) the +L-stable normal relaxation. Lowers throughflow + improves mesh quality. + +### 3. Free-slip inner boundary — rotation nullspace (`--inner freeslip`) +The rigid rotation `[-y,x]` is a velocity nullspace ONLY while the boundary is CIRCULAR. +Once the free surface DEFORMS, do NOT attach it to the held/consistent solves (→ held +22 s/`DIVERGED_LINEAR_SOLVE`, throughflow blow-up). Keep `petsc_use_pressure_nullspace`; +strip the gauge with the exact post-solve projection `_project_out_rotation` on `v`, +`v_cons` (drives advection) AND `v_h` (one consistent non-rotating frame). The undeformed +free-slip *reference* (`--surface freeslip`) is fine WITH the nullspace attached. + +### 4. Graded / higher-order meshes (drive node movement consistently) +- **Surface-ring detection**: tie the tolerance to the FINEST cell + (`0.5·mesh.get_min_radius()`), NOT the nominal `cellsize`. On a gmsh-graded mesh + (`cellSizeOuter`) the old tolerance scoops the first interior ring → a 2%-thick + "surface band" → tangling (looks like the surface "destroying itself"; it isn't — + the diffuser was fed a corrupt surface). BETTER (TODO): build the ring from the DMPlex + `Upper` label (`dm.getLabel("Upper").getStratumIS`), removing the tolerance entirely. +- **Node movement**: the solve velocity is P2; the mesh geometry is P1. Drive `u_n` and + the tangential transport from a P1 length-smoothed `Vector_Projection` of V (`v_p1`), + NOT a point-evaluation of the P2 field. +- **Stress smoothing**: `topo_proj.smoothing_length` = a fixed PHYSICAL length + (`--smooth-length`), not cell-count, so `h_∞` is mesh/order-independent. + +### 5. Cost — there is no acceleration win (don't chase it) +3 Stokes solves/step (free→u_n, held→h_∞, consistent→advect). Warm-start does NOT help +(outer KSP already 1 iter; FMG supplies its own nested guess — measured SLOWER). Blend- +skip rarely fires (α-spread always large). Operator/PC reuse: already reused across +`solve()`s (first 5.5 s setup, steady ~945 ms = irreducible FMG solve; RHS-only resolve +same cost). The per-step cost is the geometric FMG hierarchy REBUILD on the deforming +mesh — intrinsic to moving meshes, "live with it." The UNIFIED-PENALTY single solver +(`penalty·(v·n̂ − V₁·n̂)·n̂`; penalty=0→free, V₁=0→held, V₁=ũ_n→consistent — held & +consistent share the matrix) is the cleanest formulation (no recompile on the constant) +but doesn't cut the irreducible solve. + +### Elastic-plate flexure `h_∞` — IMPLEMENTED (`--flexure-D`) +Generalizes the LOCAL Airy `h_∞ = −σ_nn/ρg` (the D=0 limit) to a flexed plate +`(D ∂_s^4 + ρg) h_∞ = −σ_nn`, solved SPECTRALLY on the ring (serial Fourier — the +feasible substitute for UW3's blocked 1D-manifold FE solve): per mode +`h_∞(m) = −σ_nn(m)/(ρg + D(m/r_o)⁴)`. `D` sets the flexural wavelength `(D/ρg)^{1/4}` +and damps short-wavelength loads — the physically-grounded, mesh-independent length- +smoothing. In `_h_inf` (h_∞-ONLY — the stable form). **PROTOTYPE — amplitude response correct +(stiffer plate → less deflection) but it does NOT low-pass the SURFACE**: filtering `h_∞` only +sets a smooth set-point; the surface still picks up short-wavelength content from the SL +tangential transport + partial relaxation. "Filter every surface number (h, ḣ, h_∞)" was +TRIED and REJECTED — filtering the GEOMETRY `h` injects a spurious smooth-the-mesh motion into +`ũ_n` → flow runs away (Vrms 50→345); filtering `ḣ` alone is stable but elevates Vrms with no +benefit. So making flexure a TRUE surface low-pass without destabilising is OPEN/hard +(`_flex_filter` helper is in place). Examples: `stagnant_lid_mode1_study/{figures/flexure_*.png, +runs/flexure_D*}`. + +### Open / next (not yet done) +- **Label-based surface ring** (replace the radial heuristic with the `Upper` stratum). +- **Flexure D calibration** to a realistic lithospheric flexural wavelength. + +### Diagnostic tools (`~/+Simulations/FreeSurface/convection/stagnant_lid_mode1_study/scripts/`) +- `heldlid_hinf_check.py` — verify `h_∞` via 4 independent free-slip enforcements × Δη sweep +- `stitch_compare.py` — side-by-side montage of per-run dirs (`--dirs a,b,c`) +- `unified_penalty_solver.py` — the one-solver penalty formulation probe +- `resolve_timing_probe.py` — repeated-solve / lag-Jacobian / reuse-PC timing + +## Failure modes — symptom → cause + +| Symptom | Cause | +|---|---| +| Surface deforms but `u_n` RUNS AWAY (e.g. u_n 42→125→285→445), cold lid leaks in, plumes punch through | **RESOLVED**: material-surface advection inconsistency — T advected with stress-free `u_n` while surface moves by relaxed `ũ_n`. Fix = `--advect-velocity consistent` (Hardening §1). NOT an h_∞ bug, NOT fixed by FSSA. | +| `held` solve 22 s / `DIVERGED_LINEAR_SOLVE`, throughflow blows up, with `--inner freeslip` | rigid-rotation `[-y,x]` attached as a nullspace on the DEFORMED (non-circular) surface — invalid. Don't attach it on the moving surface; use the post-solve projection (Hardening §3). | +| Graded-mesh surface "destroys itself" (q→0.2, h_max 2% at step 1) | surface-detection tolerance scooped the first interior ring → 2%-thick band, NOT real deformation. Tie tolerance to finest cell (Hardening §4). The diffuser is innocent. | +| Stress-free top but surface not updated each step | nothing stops throughflow (the stress-free top is an open boundary unless the integrator moves the surface to track `u_n`) | +| Nu decays when it should be steady (kinematic free surface) | LAG: the SL foot reaches beyond an under-moved surface → cold pump. Fix = the h_∞ relaxation, not more smoothing | +| Surface "mountain" / one-step spike on adaptive mesh | held-lid Nitsche penalty over-stiffened by GLOBAL h; use `local_h=True` (default, PR #275) = `mesh.cell_size()` | + +## Dead ends (already tried — do NOT repeat) + +- **FSSA** signed-traction free-surface: diverges / under-deforms — rejected. +- **High-k post-smoothing** of the surface: the instability is low-m, smoothing + the wrong band. +- Per-node freeze-clamp in the relaxation (the old `relax` fatal bug). + +## Diagnose by + +`h_max` (deflection, as % of r_o), `u_n` / `vhmax` (surface throughflow — should NOT +grow unbounded), `hinf_max` (the equilibrium target), `vrms`, `Nu`. Compare the +free-surface run against the free-slip RIGID-top reference at matched physical time. diff --git a/.claude/skills/plasticity-solvers/SKILL.md b/.claude/skills/plasticity-solvers/SKILL.md new file mode 100644 index 00000000..9c0c243b --- /dev/null +++ b/.claude/skills/plasticity-solvers/SKILL.md @@ -0,0 +1,131 @@ +--- +name: plasticity-solvers +description: How to get hard-Min viscoplastic / visco-elastic-plastic (VEP) Stokes solves to CONVERGE in Underworld3 — the consistent Newton tangent (solver.consistent_jacobian), the δ-soft-min yield law, and the yield homotopy that pairs them. Reach for THIS first when a Drucker-Prager / yield-stress Stokes solve stalls, diverges (DIVERGED_LINEAR_SOLVE / line-search fail), or grinds through ~20+ nonlinear iterations. Tells you which tangent to use per model, how to confirm you are actually running Newton, and the measured failure modes. +--- + +# plasticity-solvers + +The workable recipe for **nonlinear convergence of yielding (viscoplastic / VEP) +Stokes** in Underworld3. Hard-`Min` yield laws have a non-differentiable kink that +breaks naive solvers; this encodes the combination that converges. + +**One call does it:** + +```python +stokes.constitutive_model = cm # ViscoPlastic / ViscoElasticPlastic / TI-VEP +cm.Parameters.yield_stress = tau_y # finite -> plasticity active +cm.enable_yield_homotopy() # <- the strategy: δ-ramp + the right tangent +stokes.solve(zero_init_guess=False) +``` + +`enable_yield_homotopy()` is the recommended default for any hard-Min solve. It picks +the Jacobian tangent per model (below) and ramps the soft-min δ→0 within one solve. + +--- + +## The two ingredients + +Yielding viscoplasticity is `η_eff = Min(η_visc, η_yield)`, +`η_yield = τ_y/(2·ε̇_II)`. The `Min` kink is what makes it hard. + +1. **Consistent Newton tangent** — `solver.consistent_jacobian`: + - `False` (default) — frozen-viscosity **Picard** tangent: contractive, globally + stable, **linear** (slow). Bit-identical to long-standing UW3. + - `True` — full **Newton** tangent (the assembled `dF1/dL` carries `∂η/∂(grad v)`): + **quadratic** near the solution; the kink can break it far from it. + - `"continuation"` — Picard→Newton α-blend (α a `constants[]` atom ramped 0→1). + +2. **δ-soft-min yield law + homotopy** — `g = 1 + ½(f-1+√((f-1)²+δ²)) − offset`, + `η_eff = η_ve/g`, `f = η_ve/η_pl`: + - **δ = 0 (default) ≡ exact `Min`** to machine precision. + - δ and the onset offset are `constants[]` atoms → δ is **runtime-rampable with no + JIT recompile**. `enable_yield_homotopy()` ramps δ from `delta_start`→0 within one + solve (residual-paced absolute schedule, via the `SNESSetUpdate` hook). The smooth + (δ>0) problem warm-starts the sharper one; δ ends at 0 so the **converged answer is + on the exact yield surface**. + +This is **problem-space** continuation (ramp the residual smooth→sharp). A smooth +Jacobian on a sharp `Min` residual is the consistent tangent of a *different* (harmonic) +problem and diverges worse than Picard — don't do that. + +--- + +## Confirm you are actually running Newton + +A consistent-Newton solve should converge **quadratically** — the residual roughly +squares each iteration and reaches ~1e-12 in 3–6 nonlinear steps. A **linear** tail +(residual dropping by a roughly constant factor over ~15–25 steps) means you are on the +Picard tangent — check `solver.consistent_jacobian is True` and that the viscosity is a +field of the unknowns, not a constant. + +Direct symbolic check that the Newton term is in the Jacobian (`dF1/dL` differs between +the frozen and unwrapped flux by exactly the `∂η/∂(grad v)` term): + +```python +import sympy +from underworld3.function.expressions import unwrap_expression +F1 = sympy.Array(stokes.F1.sym) # residual flux (η wrapped) +L = sympy.Array(stokes.Unknowns.L) # velocity-gradient symbols +G_picard = sympy.derive_by_array(F1, L) +F1_unwrapped = sympy.Array( + [unwrap_expression(e, mode="symbolic_keep_constants") for e in F1], F1.shape) +G_newton = sympy.derive_by_array(F1_unwrapped, L) +assert sympy.simplify(sympy.Array(G_newton) - sympy.Array(G_picard)) != \ + sympy.Array([0]*len(list(sympy.flatten(F1)))*len(list(sympy.flatten(L)))) +# nonzero difference == the Newton form is present +``` + +--- + +## Which tangent for which model (measured) + +`enable_yield_homotopy(consistent_tangent="auto")` is **model-aware**: + +| Model | `"auto"` picks | Why | +|-------|----------------|-----| +| `ViscoPlasticFlowModel` (non-elastic) | **Newton** + δ-ramp | δ-ramp keeps the residual smooth while Newton finds the basin; both sharpen to exact Min together. From a viscous guess the consistent tangent is robust (≈3–4 quadratic iters). | +| `ViscoElasticPlasticFlowModel` (VEP) | **Picard** + δ-ramp | The consistent yield tangent over the elastic stress-history block makes the Jacobian **indefinite → `DIVERGED_LINEAR_SOLVE`**. Picard is contractive; with the δ-ramp it converges to the exact yield surface. | +| `TransverseIsotropicVEPFlowModel` (TI-VEP) | **Picard** + δ-ramp | same as VEP (elastic). | + +Override with `consistent_tangent=True` / `False` / `"continuation"`. + +> Measured: VEP loading-through-yield — Picard+δ-ramp converges (σ locks at τ_y), +> Newton+δ-ramp diverges every step (`DIVERGED_LINEAR_SOLVE`). + +--- + +## Failure modes → fixes + +| Symptom | Cause | Fix | +|---------|-------|-----| +| `DIVERGED_LINEAR_SOLVE`, 0 iters, VEP | consistent Newton tangent on the elastic block → indefinite | use Picard (`consistent_tangent=False`, or `"auto"` on a VEP model) | +| Line-search failure at first step | `bt` line search trips on the kink as δ steps down | `enable_yield_homotopy` sets `snes_linesearch_type="basic"` (full step); keep it | +| Converges but σ sits **below** τ_y | solving a *smoothed* problem (fixed δ>0, or a smooth surrogate) | use the δ-ramp (δ→0) so the endpoint is exact Min | +| Linear (~20-iter) convergence | running Picard when you wanted Newton | `consistent_jacobian=True` on a non-elastic model (see "Confirm" above) | + +--- + +## Gotchas + +- **`./uw build` → `amr-dev` env.** Verify `uw.__file__` is the worktree site-packages. +- **Run VEP tests UNFORKED** — `pytest --forked` SIGABRTs here (fork of multithreaded PETSc). +- `enable_yield_homotopy` raises the `snes_max_it` floor to 100 and needs the model + **attached to a solver first** (`solver.constitutive_model = model`). +- `harmonic` yield mode is a **distinct physical model** (parallel blend), not an + approximation to Min — the homotopy does not apply to it. +- If you project η, use a **low-order** field (P0/P1) — higher order overshoots and η + is not guaranteed positive. + +--- + +## Reference + +- `ViscousFlowModel._combine_yield`, `enable_yield_homotopy` / `_yield_homotopy_step` + in `constitutive_models.py`; `solver.consistent_jacobian` / `_jacobian_source`. +- Design: `docs/developer/design/jacobian-consistent-tangent.md`. +- Tests: `tests/test_1053_yield_homotopy.py`. +- Benchmark: `docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py`. + +Footnote: before this work UW3 differentiated the flux with the viscosity still +wrapped, so `∂η/∂(grad v)` was dropped and viscoplastic solves silently ran the Picard +tangent — the origin of the "~20 iterations is intrinsic" folklore. diff --git a/.claude/skills/uw-visualisation/SKILL.md b/.claude/skills/uw-visualisation/SKILL.md new file mode 100644 index 00000000..87cba71f --- /dev/null +++ b/.claude/skills/uw-visualisation/SKILL.md @@ -0,0 +1,116 @@ +--- +name: uw-visualisation +description: Render Underworld3 mesh fields (T, V, viscosity, the adapted mesh) correctly with PyVista. Use whenever you need to SEE a UW3 result — a field colormap, the moving/adapted mesh, streamlines, or compare runs. Reach for THIS before hand-rolling a renderer; getting the four cosmetic settings wrong makes renders look grey/patchy/blocky and wastes a round-trip with Louis. +--- + +# uw-visualisation + +Canonical PyVista recipe for Underworld3 fields. This exists because every fresh +Claude session re-derives the renderer and gets the colormap / background / +lighting / DOF-sampling wrong, producing "grey/patchy/weird" images Louis +rejects. The settings below match his reference renders exactly. + +**Use PyVista (`underworld3.visualisation`), NOT matplotlib.** Louis reaffirmed +this even after seeing the legacy matplotlib renderer +(`scripts/fault_convection_frames.py`) — that one is NOT preferred. + +## Hard rules (artifacts + output location) + +- **Outputs go under `~/+Simulations/...`, NEVER `/tmp`** (Louis can't view /tmp + or harness task paths). Mirror the run's `--sim-dir`; write `T_.png` into + the run directory, comparison figures into the sim-dir root. +- `pv.OFF_SCREEN = True` at import; finish with `pl.screenshot(path); pl.close()`. + +## The field+mesh pattern (copy this exactly) + +```python +import numpy as np, underworld3 as uw, underworld3.visualisation as vis, pyvista as pv +pv.OFF_SCREEN = True + +mesh = uw.discretisation.Mesh(f"{label}.mesh.00000.h5") # or the live mesh +T = uw.discretisation.MeshVariable("T_v2p1", mesh, 1, degree=3, continuous=True) +T.read_timestep(label, "T_v2p1", 0, outputPath=D) # or use the live var + +pv_T = vis.meshVariable_to_pv_mesh_object(T) # Delaunay through T's OWN DOFs +pv_T.point_data["T"] = np.asarray(T.data[:, 0]) # attach DOF values DIRECTLY (P3-faithful) +edges = vis.mesh_to_pv_mesh(mesh).extract_all_edges() + +pl = pv.Plotter(off_screen=True, window_size=(1000, 1000)) +pl.set_background("white") # rule 2 +pl.add_mesh(pv_T, scalars="T", cmap="RdBu_r", clim=(0, 1), # rules 1 + clim required + show_edges=False, lighting=False) # rule 3 +pl.add_mesh(edges, color="black", line_width=0.5, lighting=False) # mesh overlay +pl.view_xy(); pl.camera.zoom(1.3) +pl.screenshot(out); pl.close() +``` + +## The four things that make renders look bad (all COSMETIC) + +1. `cmap="coolwarm"` → muddy grey-lavender midtone — this IS the "blue/grey/red" + Louis rejects. **Use `cmap="RdBu_r"`** (clean blue→white→red). +2. PyVista's default grey background bleeds through RdBu_r's white (T≈0.5) → dirty + grey. **Always `pl.set_background("white")`.** +3. Default lighting darkens the colormap. **Always `lighting=False`** on every + `add_mesh`. +4. Re-evaluating via `scalar_fn_to_pv_points` / vertex-only sampling drops the + high-order DOFs → blocky. **Attach `T.data[:,0]` directly** to the DOF-cloud + mesh from `meshVariable_to_pv_mesh_object` (it is correct for annulus/box/disc + — do NOT avoid it). `clim` MUST be passed (default `clim=""` trips `np.any`). + +## Seeing the MESH (adaptation / moving mesh) + +The full-annulus T colormap **washes out mesh detail** — at whole-domain zoom the +grading is invisible. To judge adaptation you MUST crop: + +- Zoom the feature region with a parallel camera: + `pl.camera.parallel_projection = True; pl.camera.parallel_scale = half_width; + pl.camera.focal_point = (cx, cy, 0)`. +- For mesh-only views, drop the field and draw `edges` on white, `line_width≈0.7`. +- Real corruption vs render artifact: apparent "holes / lumps" are often a + mesh-overlay/low-res artifact. Before calling adaptation broken, CHECK the + field's value range is bounded and count folded elements (negative cell area) + programmatically — do NOT diagnose from a render alone. +- **Overlay the feature you're refining to** (a fault trace, an interface): draw it + as a red `pv.PolyData` line over the mesh. Without it you cannot tell whether the + refinement sits ON the feature or has drifted off it (a real failure mode — see + the `adaptive-meshing` skill). Read the geometry from the run manifest so any run + renders the same way. + +## Adaptive / long runs + +- **Render each checkpoint as it lands**, not just the last frame: arm a Monitor that + polls for new `run.mesh.NNNNN.{xdmf,h5}` and emits the index → render on each + event. A completion-only watch leaves you blind for a multi-hour (e.g. TI) run. +- The per-step mesh GEOMETRY must have been written (`write_timestep(..., + meshUpdates=True)`) or you'll render deformed fields on the stale step-0 mesh. + Load the per-step `run.mesh.NNNNN.h5` as the mesh, then `read_timestep` the vars. + +## Velocity + +Same pattern; use **streamlines, not glyphs**. Build a pv mesh for V, add +`pv_mesh.streamlines(...)` or evaluate V on a line seed. Magnitude with the same +white-bg / lighting=False rules. + +## Quantities to judge a convection run (not just pretty pictures) + +- `vrms` from `uw.function.evaluate(V.sym.dot(V.sym), mesh.X.coords)` → the clean + kinetic-energy indicator (more reliable than nodal boundary metrics). +- Surface heat flux Nu via `uw.maths.BdIntegral` on the Upper boundary. +- Mesh quality: fault/bulk nearest-neighbour spacing RATIO (cKDTree) for refinement; + folded-element count + min cell area for tangling. + +## Templates in this skill + +- `render_field.py` — single/`--all`-steps T+mesh render of a run directory. +- `render_field_streamlines.py` — T colormap + mesh + **V streamlines** (sparse + seeds, thin lines, short integration so weak/closed cells read clearly, not + black spiral-blobs). Use for convection. `--tag --all`. +- `zoom_compare.py` — side-by-side cropped mesh+field for N runs at one step. + +Copy these into the run's `scripts/` (or run in place), point `--sim-dir` at the +run, and adjust the field/variable names. They already encode every rule above. + +## Related memory + +`feedback_use_uw_pyvista_visualisation.md`, `feedback_pyvista_viz_pattern.md`, +`feedback_render_all_steps.md`, `project_adaptation_corruption_was_render_artifact.md`. diff --git a/.claude/skills/uw-visualisation/render_field.py b/.claude/skills/uw-visualisation/render_field.py new file mode 100644 index 00000000..c02e29c8 --- /dev/null +++ b/.claude/skills/uw-visualisation/render_field.py @@ -0,0 +1,57 @@ +"""Render a UW3 scalar MeshVariable + mesh edges for one or all snapshots of a +run directory, following the canonical UW3 PyVista pattern (see SKILL.md). + +Outputs _.png into the run directory. Outputs live under ~/+Simulations. + +Usage: + python render_field.py --tag --sim-dir ~/+Simulations/ [--step stepNNNN | --all] + python render_field.py --tag myrun --var T_v2p1 --degree 3 --clim 0 1 --all +""" +import os, glob, re, argparse +import numpy as np, underworld3 as uw, underworld3.visualisation as vis, pyvista as pv + +pv.OFF_SCREEN = True +ap = argparse.ArgumentParser() +ap.add_argument("--tag", required=True) +ap.add_argument("--sim-dir", default="~/+Simulations/StagnantLid") +ap.add_argument("--var", default="T_v2p1", help="MeshVariable name in the checkpoint") +ap.add_argument("--degree", type=int, default=3) +ap.add_argument("--clim", type=float, nargs=2, default=[0.0, 1.0]) +ap.add_argument("--cmap", default="RdBu_r") +ap.add_argument("--step", default="", help="stepNNNN; empty = latest") +ap.add_argument("--all", action="store_true", help="render every snapshot") +args = ap.parse_args() +D = os.path.expanduser(os.path.join(args.sim_dir, args.tag)) + +cands = sorted(glob.glob(os.path.join(D, "step*.mesh.00000.h5")), + key=lambda c: int(re.search(r"step(\d+)\.mesh", c).group(1))) +if args.all: + labels = [re.search(r"(step\d+)\.mesh", os.path.basename(c)).group(1) for c in cands] +elif args.step: + labels = [args.step] +else: + labels = [re.search(r"(step\d+)\.mesh", os.path.basename(cands[-1])).group(1)] + + +def render(label): + mesh = uw.discretisation.Mesh(os.path.join(D, f"{label}.mesh.00000.h5")) + var = uw.discretisation.MeshVariable(args.var, mesh, 1, degree=args.degree, + continuous=True) + var.read_timestep(label, args.var, 0, outputPath=D) + pv_v = vis.meshVariable_to_pv_mesh_object(var) + pv_v.point_data["f"] = np.asarray(var.data[:, 0]) + edges = vis.mesh_to_pv_mesh(mesh).extract_all_edges() + pl = pv.Plotter(off_screen=True, window_size=(1000, 1000)) + pl.set_background("white") + pl.add_mesh(pv_v, scalars="f", cmap=args.cmap, clim=tuple(args.clim), + show_edges=False, lighting=False, + scalar_bar_args={"title": args.var}) + pl.add_mesh(edges, color="black", line_width=0.5, lighting=False) + pl.view_xy(); pl.camera.zoom(1.3) + out = os.path.join(D, f"{args.var}_{label}.png") + pl.screenshot(out); pl.close() + print("->", out, flush=True) + + +for lab in labels: + render(lab) diff --git a/.claude/skills/uw-visualisation/render_field_streamlines.py b/.claude/skills/uw-visualisation/render_field_streamlines.py new file mode 100644 index 00000000..c0fc4752 --- /dev/null +++ b/.claude/skills/uw-visualisation/render_field_streamlines.py @@ -0,0 +1,80 @@ +"""Render a UW3 run: T colormap + mesh edges + VELOCITY STREAMLINES. +Follows the canonical UW3 PyVista pattern (RdBu_r, white bg, lighting=False, +DOF-faithful) and overlays streamlines of V. Outputs _.png. + +Usage: + python render_field_streamlines.py --tag --sim-dir ~/+Simulations/ [--step stepNNNN | --all] +""" +import os, glob, re, argparse +import numpy as np, underworld3 as uw, underworld3.visualisation as vis, pyvista as pv +pv.OFF_SCREEN = True + +ap = argparse.ArgumentParser() +ap.add_argument("--tag", required=True) +ap.add_argument("--sim-dir", default="~/+Simulations/StagnantLid") +ap.add_argument("--tvar", default="T_v2p1"); ap.add_argument("--tdeg", type=int, default=3) +ap.add_argument("--vvar", default="V_v2p1"); ap.add_argument("--vdeg", type=int, default=2) +ap.add_argument("--clim", type=float, nargs=2, default=[0.0, 1.0]) +ap.add_argument("--step", default=""); ap.add_argument("--all", action="store_true") +ap.add_argument("--rin", type=float, default=0.5); ap.add_argument("--rout", type=float, default=1.0) +args = ap.parse_args() +D = os.path.expanduser(os.path.join(args.sim_dir, args.tag)) + +cands = sorted(glob.glob(os.path.join(D, "step*.mesh.00000.h5")), + key=lambda c: int(re.search(r"step(\d+)\.mesh", c).group(1))) +if args.all: + labels = [re.search(r"(step\d+)\.mesh", os.path.basename(c)).group(1) for c in cands] +elif args.step: + labels = [args.step] +else: + labels = [re.search(r"(step\d+)\.mesh", os.path.basename(cands[-1])).group(1)] + + +def render(label): + mesh = uw.discretisation.Mesh(os.path.join(D, f"{label}.mesh.00000.h5")) + T = uw.discretisation.MeshVariable(args.tvar, mesh, 1, degree=args.tdeg, continuous=True) + T.read_timestep(label, args.tvar, 0, outputPath=D) + V = uw.discretisation.MeshVariable(args.vvar, mesh, mesh.dim, degree=args.vdeg, continuous=True) + V.read_timestep(label, args.vvar, 0, outputPath=D) + + pv_T = vis.meshVariable_to_pv_mesh_object(T) + pv_T.point_data["T"] = np.asarray(T.data[:, 0]) + edges = vis.mesh_to_pv_mesh(mesh).extract_all_edges() + + # velocity mesh + 3D vectors for streamlines + pv_V = vis.meshVariable_to_pv_mesh_object(V) + Vd = np.asarray(V.data) + vec = np.zeros((pv_V.n_points, 3)); vec[:, 0] = Vd[:, 0]; vec[:, 1] = Vd[:, 1] + pv_V["V"] = vec; pv_V.set_active_vectors("V") + vmax = float(np.linalg.norm(Vd, axis=1).max()) + + # seed a sparse disc of points across the annulus (sparse + short + # integration so streamlines read as flow direction, not wound-up spirals + # that black out the field — low-Ra cells are tight closed loops). + rng = np.linspace(args.rin + 0.05, args.rout - 0.05, 4) + th = np.linspace(0, 2 * np.pi, 14, endpoint=False) + R, TH = np.meshgrid(rng, th) + seed = pv.PolyData(np.c_[(R * np.cos(TH)).ravel(), + (R * np.sin(TH)).ravel(), + np.zeros(R.size)]) + strm = pv_V.streamlines_from_source( + seed, vectors="V", integration_direction="both", + max_step_length=0.5, initial_step_length=0.02, + max_steps=300, terminal_speed=max(vmax * 1e-4, 1e-9)) + + pl = pv.Plotter(off_screen=True, window_size=(1000, 1000)) + pl.set_background("white") + pl.add_mesh(pv_T, scalars="T", cmap="RdBu_r", clim=tuple(args.clim), + show_edges=False, lighting=False, scalar_bar_args={"title": "T"}) + pl.add_mesh(edges, color="grey", line_width=0.3, lighting=False, opacity=0.4) + if strm.n_points > 0: + pl.add_mesh(strm, color="black", line_width=1.4, lighting=False) + pl.add_text(f"{label} |v|max={vmax:.2f}", font_size=10, color="black") + pl.view_xy(); pl.camera.zoom(1.3) + out = os.path.join(D, f"TV_{label}.png") + pl.screenshot(out); pl.close() + print("->", out, flush=True) + + +for lab in labels: + render(lab) diff --git a/.claude/skills/uw-visualisation/zoom_compare.py b/.claude/skills/uw-visualisation/zoom_compare.py new file mode 100644 index 00000000..0025aaab --- /dev/null +++ b/.claude/skills/uw-visualisation/zoom_compare.py @@ -0,0 +1,59 @@ +"""Zoomed side-by-side of the fault region (12 o'clock) showing T + the +adapted mesh for several runs at one step. The full-annulus render washes out +fault mesh detail (memory note); this crops to the fault neighbourhood so the +clustering (or lack of it) is visible. + +Usage: + python fault_zoom_compare.py --step step0060 \ + --tags rq_passive_uniform rq_passive_gmsh --labels "uniform 1.6x" "gmsh 4.7x" +""" +import os, argparse +import numpy as np, underworld3 as uw, underworld3.visualisation as vis, pyvista as pv + +pv.OFF_SCREEN = True +ap = argparse.ArgumentParser() +ap.add_argument("--step", default="step0060") +ap.add_argument("--tags", nargs="+", required=True) +ap.add_argument("--labels", nargs="+", default=None) +ap.add_argument("--sim-dir", default="~/+Simulations/StagnantLid+Fault") +ap.add_argument("--out", default="fault_zoom_compare.png") +# Fault is at theta=90 (top), dipping east. Crop a box around the trace. +ap.add_argument("--cx", type=float, default=0.12) +ap.add_argument("--cy", type=float, default=0.82) +ap.add_argument("--half", type=float, default=0.34) +args = ap.parse_args() + +SIM = os.path.expanduser(args.sim_dir) +labels = args.labels or args.tags +n = len(args.tags) + +pl = pv.Plotter(off_screen=True, shape=(1, n), window_size=(700 * n, 760), + border=False) +pl.set_background("white") +for i, (tag, lab) in enumerate(zip(args.tags, labels)): + D = os.path.join(SIM, tag) + mp = os.path.join(D, f"{args.step}.mesh.00000.h5") + if not os.path.exists(mp): + continue + mesh = uw.discretisation.Mesh(mp) + T = uw.discretisation.MeshVariable("T_v2p1", mesh, 1, degree=3, + continuous=True, varsymbol="T") + T.read_timestep(args.step, "T_v2p1", 0, outputPath=D) + pv_T = vis.meshVariable_to_pv_mesh_object(T) + pv_T.point_data["T"] = np.asarray(T.data[:, 0]) + edges = vis.mesh_to_pv_mesh(mesh).extract_all_edges() + pl.subplot(0, i) + pl.add_text(lab, font_size=11, color="black") + pl.add_mesh(pv_T, scalars="T", cmap="RdBu_r", clim=(0, 1), + show_edges=False, lighting=False, show_scalar_bar=False) + pl.add_mesh(edges, color="black", line_width=0.7, lighting=False) + pl.view_xy() + pl.camera.focal_point = (args.cx, args.cy, 0.0) + pl.camera.position = (args.cx, args.cy, 10.0) + pl.camera.parallel_projection = True + pl.camera.parallel_scale = args.half + +out = os.path.join(SIM, args.out) +pl.screenshot(out) +pl.close() +print("->", out) diff --git a/.gitignore b/.gitignore index 7e962177..18fb2948 100644 --- a/.gitignore +++ b/.gitignore @@ -204,12 +204,21 @@ Jupyterbook/Notebooks/.DS_Store Jupyterbook/Notebooks/Examples-NavierStokes/order_1/NS_test_Re_250_0.06_step_0.png Documentation/index.aux -# Claude Code files (track commands, ignore settings) +# Claude Code files (track commands and skills, ignore settings) CLAUDE.md .claude/* !.claude/commands/ .claude/commands/* !.claude/commands/*.md +# Ship skills with the repo: allow SKILL.md + supporting assets, ignore junk +!.claude/skills/ +.claude/skills/**/* +!.claude/skills/**/ +!.claude/skills/**/*.md +!.claude/skills/**/*.py +!.claude/skills/**/*.typ +!.claude/skills/**/*.json +!.claude/skills/**/*.png # Prevent future markdown clutter in root (allow essential docs only) /*_ANALYSIS.md