Simplify Priority 1 APIs around functions

AGENTS.md

@@ -47,6 +47,7 @@ geodef/
 | `greens` | Green's matrix assembly, projection, stacking, Laplacian operators |
 | `gradients` | Differentiable forward models: Jacobians w.r.t. geometry and slip (JAX) |
 | `fault` | `Fault` class: factory methods, forward modeling, I/O, moment |
+| `slip` | Slip-vector packing and strike/dip, rake, azimuth, and plate-basis conversions |
 | `medium` | `ElasticMedium`: shear modulus and Poisson's ratio, shared by Green's functions, stress kernels, and moment |
 | `data` | `DataSet` base + `GNSS`, `InSAR`, `Vertical` data types |
 | `invert` | Inversion: solvers, fixed-direction slip bases, regularization, hyperparameter tuning, model assessment, scalar/per-component/per-parameter bounds |

CLAUDE.md

@@ -47,6 +47,7 @@ geodef/
 | `greens` | Green's matrix assembly, projection, stacking, Laplacian operators |
 | `gradients` | Differentiable forward models: Jacobians w.r.t. geometry and slip (JAX) |
 | `fault` | `Fault` class: factory methods, forward modeling, I/O, moment |
+| `slip` | Slip-vector packing and strike/dip, rake, azimuth, and plate-basis conversions |
 | `medium` | `ElasticMedium`: shear modulus and Poisson's ratio, shared by Green's functions, stress kernels, and moment |
 | `data` | `DataSet` base + `GNSS`, `InSAR`, `Vertical` data types |
 | `invert` | Inversion: solvers, fixed-direction slip bases, regularization, hyperparameter tuning, model assessment, scalar/per-component/per-parameter bounds |

README.md

@@ -85,18 +85,18 @@ ue, un, uz = fault.displacement(obs_lat, obs_lon, slip_strike=0.0, slip_dip=1.0)
 gnss = GNSS.load("stations.dat")
 insar = geodef.InSAR.load("ascending.dat")
 
-result = geodef.invert(fault, [gnss, insar],
-                       smoothing='laplacian',
-                       smoothing_strength=1e3,
-                       bounds=(0, None))
+result = geodef.solve(fault, [gnss, insar],
+                      smoothing='laplacian',
+                      smoothing_strength=1e3,
+                      bounds=(0, None))
 
 print(f"Mw = {result.Mw:.2f}, reduced chi2 = {result.reduced_chi2:.2f}")
 geodef.plot.slip(fault, result.slip_vector)
 
 # Optional fixed slip directions
-fixed_rake = geodef.invert(fault, gnss, components='rake', rake=90.0)
-fixed_azimuth = geodef.invert(fault, gnss,
-                              components='azimuth', slip_azimuth=15.0)
+fixed_rake = geodef.solve(fault, gnss, components='rake', rake=90.0)
+fixed_azimuth = geodef.solve(fault, gnss,
+                             components='azimuth', slip_azimuth=15.0)
 ```
 
 ## Differentiable and Bayesian modeling (JAX)
@@ -164,6 +164,7 @@ Full API docs with examples are in `docs/`:
 | Doc | Module |
 |-----|--------|
 | [`docs/fault.md`](docs/fault.md) | `Fault` class — factory methods, forward modeling, I/O |
+| [`docs/slip.md`](docs/slip.md) | Slip-vector functions, plate-motion coordinates, and patch ordering |
 | [`docs/medium.md`](docs/medium.md) | `ElasticMedium` half-space parameters |
 | [`docs/data.md`](docs/data.md) | `GNSS`, `InSAR`, `Vertical` data types |
 | [`docs/greens.md`](docs/greens.md) | Green's matrix assembly and Laplacian operators |

docs/backend.md

@@ -100,8 +100,8 @@ and rerun the final inversion in float64:
 geodef.backend.set_precision("float32")
 ac = geodef.abic_curve(fault, data, smoothing="laplacian")   # fast sweep
 geodef.backend.set_precision("float64")
-result = geodef.invert(fault, data, smoothing="laplacian",
-                       smoothing_strength=ac.optimal)         # final solve
+result = geodef.invert.solve(fault, data, smoothing="laplacian",
+                             smoothing_strength=ac.optimal)   # final solve
 ```
 
 With the JAX backend, precision is synced to JAX's `jax_enable_x64` flag,

docs/bayes.md

@@ -46,8 +46,8 @@ models with small `L m`; for a Laplacian, that means neighboring patches tend
 to have similar slip. These assumptions are choices, not universal physical
 laws, and should be reported with the result.
 
-- `theta` — planar-fault geometry `[e0, n0, depth, strike, dip, length,
-  width]`; any subset can be sampled (`free`), the rest stay fixed.
+- `geometry` — a parameter mapping or expert/JAX vector; any subset of its
+  named parameters can be sampled (`free`), while the rest stay fixed.
 - `sigma` — dimensionless noise scale factor multiplying the dataset
   covariances (`sigma = 1` means the reported data errors are exact).
   Sampled as `log10_sigma`.
@@ -80,10 +80,17 @@ just the best-fitting slip model.
 ## Building a posterior
 
 ```python
+frame = geodef.LocalFrame(-2.0, 100.0, projection="wgs84-enu")
+geometry0 = {
+    "e0": 0.0, "n0": 0.0, "depth": 25e3,
+    "strike": 315.0, "dip": 30.0,
+    "length": 180e3, "width": 90e3,
+}
+
 post = geodef.bayes.RectPosterior(
-    theta0,                       # [e0, n0, depth, strike, dip, length, width]
+    geometry0,
     [gnss, insar],                # any DataSet mix
-    ref_lat=-2.0, ref_lon=100.0,  # local Cartesian anchor
+    frame=frame,
     free=["dip", "depth"],
     theta_prior={
         "dip": (5.0, 60.0),               # uniform
@@ -98,9 +105,14 @@ post = geodef.bayes.RectPosterior(
 
 post.param_names   # ['dip', 'depth', 'log10_sigma', 'log10_lambda']
 post.x0            # starting point (theta0 values + scale defaults)
+post.geometry(post.x0)  # named parameter mapping for one state
+post.fault(post.x0)     # concrete Fault for forward modeling or plotting
 post.logpdf(x)     # traceable, differentiable log-posterior
 ```
 
+The seven-element `theta0` array is also supported with either
+`frame=frame` or `ref_lat=..., ref_lon=...`.
+
 `logpdf`, `log_likelihood`, and `log_prior` are pure JAX-traceable
 functions of the sampled vector `x` — hand them to any JAX sampler or
 optimizer. Uniform-prior parameters are clipped to their bounds before
@@ -339,9 +351,8 @@ enough to set priors on directly.
 ### Setup workflow — look before you sample
 
 ```python
-fault = geodef.Fault.from_triangles(nodes, ref_lat, ref_lon, triangles=tri)
-
-warp = geodef.bayes.TriWarp(fault, n_knots=(3, 2))   # or knots=(nk, 2) array
+tri_fault = geodef.Fault.from_triangles(nodes, triangles=tri, frame=frame)
+warp = geodef.bayes.TriWarp(tri_fault, n_knots=(3, 2))
 warp.knots_uv, warp.knots_xyz    # where the knots sit
 warp.length_scale                # RBF smoothness (m); override if needed
 

docs/cache.md

@@ -39,7 +39,7 @@ geodef.cache.clear()    # delete all cached files
 
 ## How caching works
 
-`geodef.greens.greens()` and `Fault.stress_kernel()` automatically cache their
+`geodef.greens.matrix()` and `Fault.stress_kernel()` automatically cache their
 results. The cache key is computed from all input arrays and parameters (fault
 geometry, observation coordinates, data class, active GNSS components, look
 vectors, etc.) using SHA-256. If the key matches an existing `.npz` file, the
@@ -67,7 +67,7 @@ To bypass caching entirely for a session:
 
 ```python
 geodef.cache.disable()
-G = geodef.greens.greens(fault, gnss)  # always recomputes
+G = geodef.greens.matrix(fault, gnss)  # always recomputes
 geodef.cache.enable()
 ```
 

docs/conventions.md

@@ -15,10 +15,13 @@ mapping is given here.
   (`Fault.centers` predates this policy and stores `[lat, lon, depth]`;
   use `Fault.centers_geo` for the documented `[lon, lat, depth]` order.)
 - **Local Cartesian:** East, North, Up (ENU) in meters, right-handed, tied
-  to an explicit reference origin (`ref_lat`, `ref_lon`). Anything named
-  `*_enu` or `*_local` uses this frame; `z`/`up` is negative below the
-  surface, while `depth` is positive down. Convert with
-  `geodef.transforms`.
+  to a `LocalFrame` that records origin latitude, longitude, altitude, and
+  projection. The current projection identifier is `"wgs84-enu"` (WGS84
+  geographic → ECEF → tangent ENU); it is explicit so later projections do
+  not make saved or combined local arrays ambiguous. Anything named `*_enu`
+  or `*_local` uses its object's `.frame`; `z`/`up` is negative below the
+  surface, while `depth` is positive down. Use `LocalFrame.to_enu`,
+  `to_geographic`, and explicit `transform_enu`/`geometry.to_frame` methods.
 - Kernel-native frames (Okada's fault-aligned x/y, DC3D internals,
   triangular dislocation coordinates) never appear in public signatures;
   they are converted at the adapter layer inside `geodef.greens`.
@@ -38,6 +41,10 @@ All angles are degrees.
 - **Slip azimuth:** geographic azimuth of horizontal slip, clockwise from
   North. `components='azimuth'` converts to each patch's local rake as
   `slip_azimuth - strike_i`, so it remains meaningful on curved meshes.
+- **Plate rake:** a large-scale kinematic direction expressed in each patch's
+  local strike/dip plane. Plate coordinates are rake-parallel and
+  rake-perpendicular; unlike raw triangle-local components, they can remain a
+  smooth basis across a variable-orientation mesh.
 
 ## Units
 
@@ -53,8 +60,11 @@ All angles are degrees.
 
 - **Slip vectors** are blocked: for `N` patches and both components, the
   first `N` entries are strike-slip, the last `N` dip-slip
-  (`m[:N]`, `m[N:]`). Single-component bases (`'strike'`, `'dip'`,
-  `'rake'`, `'azimuth'`) have length `N`.
+  (`m[:N]`, `m[N:]`). Plate coordinates are likewise blocked
+  `[rake_parallel | rake_perpendicular]`. Single-component bases (`'strike'`,
+  `'dip'`, `'rake'`, `'azimuth'`) have length `N`. Use `geodef.slip` conversion
+  functions or the named arrays on `InversionResult` to recover physical
+  strike/dip components.
 - **Green's matrix rows** follow each dataset's observation vector:
   GNSS with 3 components interleaves `[E, N, U]` per station (`[E, N]` for
   2-component data); InSAR contributes one LOS row per pixel; `Vertical`
@@ -66,6 +76,8 @@ All angles are degrees.
 - **Patch order** for structured rectangular grids varies along strike
   fastest: patch `k = i_strike + n_length * j_dip`; use
   `Fault.patch_index(strike_idx, dip_idx)` instead of hand-computing this.
+  `Fault.reshape_patches` converts patch-first arrays to
+  `[dip_index, strike_index, ...]`; `Fault.flatten_patches` reverses it.
 
 ## Regularization
 

docs/data.md

@@ -158,7 +158,7 @@ insar = InSAR(lon=lon, lat=lat, los=los, sigma=sigma, look_e=look_e,
 `load()` reads diagonal-uncertainty files. If a loaded dataset needs a full
 covariance, reconstruct it from the source arrays and pass `covariance=`.
 
-The covariance is used automatically by `geodef.invert()` and
+The covariance is used automatically by `geodef.invert.solve()` and
 `geodef.stack_weights()`.
 
 ### Building a spatially-correlated covariance

docs/fault.md

@@ -21,12 +21,9 @@ Up and uses depth positive downward.
 
 ## Factory classmethods
 
-### `Fault.planar(*, lat, lon, depth, strike, dip, length, width, n_length=1, n_width=1, medium=None)`
+### `Fault.planar(*, lat, lon, depth, strike, dip, length, width, n_length=1, n_width=1, medium=None, frame=None)`
 
-All arguments are keyword-only so latitude/longitude order can never be
-confused (see [`conventions.md`](conventions.md)).
-
-Create a discretized planar fault from its centroid.
+Create a discretized planar fault directly from named geographic parameters:
 
 ```python
 fault = Fault.planar(
@@ -38,6 +35,15 @@ fault = Fault.planar(
 # → Fault with 50 rectangular patches, engine='okada'
 ```
 
+Geographic arguments are keyword-only so latitude/longitude order cannot be
+confused. Pass an explicit frame when the fault must share local coordinates
+with another object:
+
+```python
+frame = geodef.LocalFrame(-2.0, 100.0)
+fault = Fault.planar(..., frame=frame)
+```
+
 ### `Fault.planar_from_corner(lat, lon, depth, strike, dip, length, width, n_length=1, n_width=1)`
 
 Same as `planar()` but the reference point is the top-left (shallowest, along-strike start) corner instead of the centroid.
@@ -52,17 +58,19 @@ mesh = from_slab2("sum_slab2_dep.grd", bounds=(95, 106, -6, 6))
 fault = Fault.from_mesh(mesh)
 ```
 
-### `Fault.from_triangles(vertices, ref_lat=0.0, ref_lon=0.0, *, triangles=None)`
+### `Fault.from_triangles(vertices, *, frame=None, ref_lat=None, ref_lon=None, triangles=None)`
 
-Create a triangular fault directly from ENU vertex coordinates. Two forms:
+Create a triangular fault from ENU arrays in either of two forms:
 
 - Explicit corners: `vertices` has shape `(N, 3, 3)` (leave `triangles=None`).
 - Node array + connectivity: pass a shared `(M, 3)` node array as `vertices`
   plus an `(N, 3)` index array as `triangles`. This preserves the exact patch
   order and node sharing of an imported mesh.
 
 ```python
-fault = Fault.from_triangles(nodes, ref_lat, ref_lon, triangles=tris)
+fault = Fault.from_triangles(
+    nodes, frame=frame, triangles=tris
+)
 ```
 
 ### `Fault.load(fname, *, format=None, ref_lat=0.0, ref_lon=0.0)`
@@ -90,6 +98,7 @@ fault = Fault.load("cascadia", format="ned")  # reads cascadia.ned + cascadia.tr
 |----------|-------|-------------|
 | `n_patches` | scalar | Number of patches |
 | `engine` | `str` | `"okada"` or `"tri"` |
+| `frame` | `LocalFrame` | Frame defining every local-coordinate view and triangular vertex |
 | `grid_shape` | `(nL, nW)` or `None` | Structured grid dimensions |
 | `centers` | `(N, 3)` | Patch centers as `[lat, lon, depth_m]` (legacy latitude-first order) |
 | `centers_geo` | `(N, 3)` | Patch centers as `[lon, lat, depth_m]` (documented geographic order, matches `Mesh.centers_geo`) |
@@ -104,18 +113,27 @@ fault = Fault.load("cascadia", format="ned")  # reads cascadia.ned + cascadia.tr
 
 All geometry arrays are read-only after construction.
 
+Use `fault.to_frame(target_frame)` to explicitly re-express local views. For
+triangular faults it transforms every vertex while preserving its geographic
+position; incompatible frames are never silently substituted.
+
 ---
 
 ## Forward modeling
 
 ### `fault.displacement(obs_lat, obs_lon, slip_strike, slip_dip=0.0)`
 
-Compute surface displacements for a slip distribution. `slip_strike` and
-`slip_dip` may be scalars broadcast to every patch or arrays with shape `(N,)`.
+Compute surface displacements from strike-slip and dip-slip scalars or arrays.
 
 ```python
-ue, un, uz = fault.displacement(obs_lat, obs_lon, slip_strike=0.0, slip_dip=1.0)
-# ue, un, uz each have shape (n_obs,)
+strike_slip = np.zeros(fault.n_patches)
+dip_slip = np.ones(fault.n_patches)
+east, north, up = fault.displacement(
+    obs_lat,
+    obs_lon,
+    slip_strike=strike_slip,
+    slip_dip=dip_slip,
+)
 ```
 
 ### `fault.greens_matrix(obs_lat, obs_lon, kind="displacement", obs_depth=None)`
@@ -142,8 +160,9 @@ well the spatial slip distribution is resolved. `slip` here is slip magnitude,
 not a signed strike- or dip-slip component.
 
 ```python
-M0 = fault.moment(slip, mu=30e9)      # slip magnitude shape (N,); returns N·m
-Mw = fault.magnitude(slip, mu=30e9)  # moment magnitude
+slip_magnitude = np.hypot(strike_slip, dip_slip)
+M0 = fault.moment(slip_magnitude, mu=30e9)     # returns N·m
+Mw = fault.magnitude(slip_magnitude, mu=30e9) # moment magnitude
 
 # Module-level utilities
 from geodef import moment_to_magnitude, magnitude_to_moment
@@ -175,6 +194,9 @@ friction convention, and normal-stress sign convention.
 ```python
 idx = fault.patch_index(strike_idx=3, dip_idx=1)
 # Only valid for structured grids (Fault.planar or Fault.load with grid)
+
+grid = fault.reshape_patches(values)  # (N, ...) -> (n_width, n_length, ...)
+values = fault.flatten_patches(grid)  # inverse conversion
 ```
 
 ---