diff --git a/AGENTS.md b/AGENTS.md
index a2fd438..22e0366 100644
--- a/AGENTS.md
+++ b/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 |
diff --git a/CLAUDE.md b/CLAUDE.md
index 25bd62a..195c30d 100644
--- a/CLAUDE.md
+++ b/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 |
diff --git a/PLAN.md b/PLAN.md
index 42cde88..32ac43c 100644
--- a/PLAN.md
+++ b/PLAN.md
@@ -2,9 +2,9 @@
 
 GeoDef's next phase is not a choice between a small teaching library and a
 powerful research package. It should be both: a novice should be able to express
-a geophysical problem in a few named objects, while an expert can still reach
-the Green's matrices, covariance operators, autodiff kernels, constrained
-solvers, and Bayesian posteriors underneath.
+a geophysical problem with a few memorable functions and durable domain values,
+while an expert can still reach the Green's matrices, covariance operators,
+autodiff kernels, constrained solvers, and Bayesian posteriors underneath.
 
 This roadmap replaces the implementation diary that previously occupied this
 file. Git history, tests, and the module documentation preserve the details of
@@ -72,6 +72,11 @@ These are acceptance criteria for every roadmap item.
 8. **Array transparency without array traps.** NumPy arrays remain accepted and
    returned where natural, but common shape/order mistakes should be prevented
    by named accessors and precise validation.
+9. **Functions are the default abstraction.** Add a public class only when it
+   owns durable state, preserves an invariant that functions cannot preserve,
+   or amortizes expensive preparation across calls. Parameter bundles and
+   transient array views stay as keyword arguments, arrays, and module
+   functions.
 
 ---
 
@@ -80,10 +85,10 @@ These are acceptance criteria for every roadmap item.
 ### What already serves learners well
 
 - The core relation `d = G m` is visible rather than hidden behind a framework.
-- `Fault.planar(...)`, `fault.displacement(...)`, and `geodef.invert(...)` form
+- `Fault.planar(...)`, `fault.displacement(...)`, and `geodef.solve(...)` form
   a compact path from geometry to a solution.
 - NumPy is the default; advanced compilation and sampling are opt-in.
-- Domain objects are immutable, synthetic tutorials are reproducible, and the
+- Durable domain objects are immutable, synthetic tutorials are reproducible, and the
   teaching sequence follows the concepts of an inverse problem rather than the
   package's module layout.
 - Low-level engines are independently accessible and extensively cross-checked,
@@ -184,7 +189,7 @@ These are acceptance criteria for every roadmap item.
 ## Priority 0 — Consistency and trust (small, high-impact work)
 
 Complete these before adding a new high-level abstraction. They establish the
-semantics that later objects will wrap.
+semantics used by the function-oriented public API.
 
 ### 0.1 Audit and freeze mathematical conventions
 
@@ -262,79 +267,88 @@ semantics that later objects will wrap.
 
 ---
 
-## Priority 1 — A coherent everyday workflow
-
-The functional API remains supported. This phase adds domain names and context;
-it does not hide the linear algebra or create a mandatory framework.
-
-### 1.1 Named geometry and coordinate frames
-
-- [ ] Add immutable `LocalFrame(origin_lat, origin_lon, origin_alt=0)` and
-  `PlanarGeometry(center, depth, strike, dip, length, width)` value objects.
-  They validate units/conventions and convert explicitly to geographic arrays,
-  ENU arrays, and the autodiff parameter vector.
-- [ ] Let `Fault.planar` and geometry-search/Bayesian constructors accept these
-  objects while preserving keyword-based scalar calls.
-- [ ] Replace unexplained seven-element `theta` arrays in beginner-facing docs
-  and results with named geometry views; keep `.theta` as the expert/JAX view.
-- [ ] Give every local-coordinate-bearing object a `.frame` and reject combining
-  objects with incompatible frames unless the user explicitly transforms them.
-
-### 1.2 One canonical slip representation
-
-- [ ] Add an immutable `SlipModel` with named per-patch fields/accessors:
-  `strike`, `dip`, `magnitude`, `rake`, and `.vector` for the blocked linear-
-  algebra view. Support one-component rake/azimuth amplitudes without pretending
-  they are already two components.
-- [ ] Accept `SlipModel` anywhere a slip vector is accepted; keep NumPy arrays
-  fully supported for low-level and backwards-compatible code.
-- [ ] Make forward results a small named `Displacement(east, north, up)` object
-  with tuple unpacking and `.vector` so existing idioms remain concise.
-- [ ] Standardize patch ordering utilities and provide `fault.reshape_patches`
+## Priority 1 — A small, function-oriented everyday API
+
+The public vocabulary should remain smaller than the set of concepts in the
+implementation. GeoDef keeps objects for durable domain state and prepared
+computations; ordinary transformations and one-shot workflows are functions.
+
+### 1.1 Set and enforce the public object budget
+
+- [x] Keep `Fault`, `GNSS`/`InSAR`/`Vertical`, `Mesh`, `LocalFrame`, and
+  `ElasticMedium`: each carries durable identity or invariants shared by many
+  computations. Keep immutable result records where named fields prevent tuple
+  or ordering mistakes.
+- [x] Remove the draft `PlanarGeometry` and `TriGeometry` public wrappers before
+  release. `Fault` already represents rectangular or triangular geometry;
+  retain the useful validation and conversions as functions.
+- [x] Replace the draft `SlipModel` and `Displacement` wrappers with `slip`
+  conversion functions, ordinary arrays, direct named `InversionResult` views,
+  and the existing three-array displacement return.
+- [x] Keep `LinearSystem` as an expert prepared/cache object for repeated
+  sweeps and assessment, not as the beginner workflow and not behind a second
+  `SlipProblem` facade.
+- [x] Require an explicit justification in API review for every new public
+  class: durable state, enforced cross-call invariant, or measured reuse of an
+  expensive preparation. Do not add classes solely to bundle keyword arguments.
+
+### 1.2 Give the functional namespaces memorable names
+
+- [x] Make the module path the primary discovery surface. Prefer specific names
+  over umbrella verbs: `geodef.invert.solve`, `lcurve`, `abic_curve`,
+  `dataset_diagnostics`, and `model_covariance`; `geodef.greens.matrix`,
+  `project`, and `laplacian`; `geodef.slip.pack`, `unpack`, `from_rake`,
+  `from_azimuth`, `from_plate`, `to_plate`, `magnitude`, and `rake`.
+- [x] Resolve the top-level `geodef.invert` function/module collision directly
+  before release: `geodef.invert` is the module, `geodef.invert.solve(...)` is
+  the primary call, and `geodef.solve(...)` is the short alias. No deprecation
+  shim or callable-module proxy is needed because the draft API has no users.
+- [x] Represent a slip basis with explicit function keywords (`components`,
+  `rake`, `slip_azimuth`, `plate_rake`) and conversion functions. Do not create
+  `SlipBasis`, `Regularization`, or `Bounds` configuration classes.
+- [x] Standardize patch ordering utilities and provide `fault.reshape_patches`
   / `fault.flatten_patches` rather than requiring learners to know which grid
   axis varies fastest.
 
-### 1.3 Turn `LinearSystem` into the reusable problem object
-
-- [ ] Prototype a beginner-named `SlipProblem` facade over `LinearSystem`, with
-  `fault`, named datasets, slip basis, regularizer, and noise models bound once.
-  Validate the name and call sequence with novice users before freezing it.
-- [ ] Provide the discoverable sequence
-  `problem.solve()`, `problem.select_regularization()`,
-  `problem.greens_matrix`, and `problem.assess(result)`. Keep
-  `geodef.invert(...)` as the shortest one-shot path and `LinearSystem` as a
-  compatible expert alias or implementation detail.
-- [ ] Replace overloaded strings and loosely related keyword groups with small
-  optional specifications (`SlipBasis`, `Regularization`, `Bounds`) while still
-  accepting today's strings/tuples.
-- [ ] Add `.describe()` / rich `repr` output summarizing data counts, parameter
-  counts, units, solver, covariance type, regularization, constraints, backend,
-  and estimated dense-memory cost before a solve.
-
-### 1.4 Results that know what they describe
-
-- [ ] Return a context-rich result that records fault identity, dataset names
-  and slices, slip basis, solver status, regularization selection, backend,
-  warnings, and version/provenance needed to reproduce the solve.
-- [ ] Add `result.prediction(dataset_or_name)`, `result.residual(dataset_or_name)`,
-  and `result.diagnostics(dataset_or_name)`; eliminate manual stacked-vector
-  slicing from all beginner documentation.
-- [ ] Add focused conveniences such as `result.plot_slip()`,
-  `result.plot_fit(dataset=...)`, and `result.summary()` by delegating to the
-  existing plotting/assessment functions, not duplicating their logic.
+### 1.3 Keep coordinates named without duplicating geometry
+
+- [x] Attach one immutable `LocalFrame` to every `Fault` and `Mesh` that owns
+  local coordinates; reject incompatible frames unless explicitly transformed.
+- [x] Put conversions in `transforms`/`geometry` functions accepting
+  `frame=...`. Keep `Fault.planar(...)` keyword geometry and
+  `Fault.from_triangles(..., frame=...)` as the named construction paths.
+- [x] Replace unexplained seven-element geometry arrays in beginner docs with
+  keyword calls and named result fields. Advanced JAX/Bayesian functions may
+  retain `theta`, but must accept a mapping keyed by `e0`, `n0`, `depth`,
+  `strike`, `dip`, `length`, and `width` in addition to the array view.
+- [x] Return optimized geometry as a `Fault` plus the expert `theta` array and
+  frame on the existing result record, rather than introducing a parallel
+  geometry value hierarchy.
+
+### 1.4 Result records plus assessment functions
+
+- [x] Keep `InversionResult` a compact, serializable data record. Add direct
+  named slip views (`strike_slip`, `dip_slip`, `slip_magnitude`, `slip_rake`)
+  where they are unambiguous; keep `slip_vector` as the blocked expert view.
+- [ ] Record dataset names and row slices, solver status, regularization
+  selection, backend, warnings, and minimal provenance needed to interpret and
+  reproduce a solve. Do not retain live `Fault` or dataset objects in results.
+- [ ] Add module functions `invert.prediction`, `invert.residual`,
+  `invert.diagnostics`, and `invert.summary`, plus corresponding `plot`
+  functions. Do not turn the result record into a workflow facade.
 - [ ] Define a versioned, safe result file schema with metadata and migration;
   retain `.npz` portability and add a human-readable manifest.
 
-### 1.5 Friendlier data ingestion
+### 1.5 Friendlier data functions
 
-- [ ] Add named constructors such as `GNSS.from_components(...)`,
-  `GNSS.horizontal(...)`, `InSAR.from_look_vector(...)`, and
-  `InSAR.from_incidence_heading(...)`; make optional vertical components truly
-  optional in the friendly path.
-- [ ] Add table ingestion with explicit column mappings, units, missing-value
-  handling, and station names. Keep dataframe libraries optional and accept the
-  Python dataframe interchange protocol rather than coupling the core to one
-  implementation.
+- [ ] Add `data.gnss`, `data.horizontal_gnss`, `data.insar`, and
+  `data.vertical` functions with keyword-only component names and sensible
+  defaults. They return the existing validated dataset classes; class
+  constructors remain available for compatibility.
+- [ ] Add `data.from_table` with explicit column mappings, units,
+  missing-value handling, and station names. Keep dataframe libraries optional
+  and accept the Python dataframe interchange protocol rather than coupling the
+  core to one implementation.
 - [ ] Separate displacement from velocity semantics in metadata (units and
   epoch/time span) without duplicating all dataset classes.
 - [ ] Introduce dataset names as first-class identifiers so joint results and
@@ -349,8 +363,8 @@ it does not hide the linear algebra or create a mandatory framework.
 - [ ] Create a five-minute, copy-paste quickstart that performs forward
   modeling, adds synthetic noise, solves slip, and plots observations versus
   predictions without manual vector packing or slicing.
-- [ ] Add a visual workflow page linking the four levels of API:
-  domain workflow → reusable problem → matrices/operators → physics kernels.
+- [ ] Add a visual workflow page linking the three levels of API:
+  domain functions → matrices/operators → physics kernels.
 - [ ] Add a glossary of geophysical and inverse-theory terms, with package names
   beside the mathematical symbols.
 - [ ] Provide “which function do I use?” and “which assumption am I making?”
@@ -437,18 +451,19 @@ not reorganize numerical reference ports merely to make their style conventional
   published sources; wrap them with clearer adapters rather than cosmetically
   rewriting formulas.
 
-### 3.3 Replace string dispatch with small protocols
+### 3.3 Replace string dispatch with callable contracts
 
-- [ ] Define minimal protocols for Green's engines, data projection, noise
-  whitening, regularization operators, and solvers. Built-ins use the same
-  protocol third parties can implement.
+- [ ] Define typed function signatures for Green's engines, data projection,
+  noise whitening, regularization operators, and solvers. Accept callables and
+  SciPy-style operators directly; introduce a public protocol type only where
+  static typing materially improves extension safety.
 - [ ] Register engines explicitly instead of expanding `if engine == ...`
   branches across `Fault`, `greens`, gradients, Bayesian code, and plotting.
 - [ ] Require engine capability declarations (surface/internal displacement,
   strain, autodiff, supported source geometry) and produce actionable errors
   when a workflow requests an unsupported capability.
 - [ ] Avoid a plugin framework until at least two external engines demonstrate
-  the protocol; start with ordinary Python objects and registration.
+  the callable contract; start with ordinary functions and registration.
 
 ### 3.4 Strengthen numerical contracts
 
@@ -469,9 +484,9 @@ not reorganize numerical reference ports merely to make their style conventional
 
 ### 4.1 Noise and whitening operators
 
-- [ ] Introduce a `NoiseModel`/whitener interface with diagonal, dense,
-  block-diagonal, sparse, low-rank-plus-diagonal, and user-supplied linear-
-  operator implementations.
+- [ ] Accept a whitening callable or `LinearOperator` and provide constructor
+  functions for diagonal, dense, block-diagonal, sparse, and low-rank-plus-
+  diagonal cases. Do not require users to adopt a `NoiseModel` class hierarchy.
 - [ ] Solve via whitening or factorizations rather than explicitly forming
   `W = C^-1`. Preserve `stack_weights()` as an educational/small-problem helper.
 - [ ] Add parametric spatial covariance fitting, variograms, and honest
@@ -529,15 +544,17 @@ manually re-signing slip vectors.
 - [ ] Make rectangular and triangular strain/stress kernels traceable,
   gradient-safe, jitted, and vmapped; validate derivatives against finite
   differences away from documented singular boundaries.
-- [ ] Remove hidden dependence on global backend state from compiled problem
-  objects while retaining `set_backend(...)` as the simple entry point.
+- [ ] Remove hidden dependence on global backend state from compiled kernels
+  and prepared systems while retaining `set_backend(...)` as the simple entry
+  point.
 - [ ] Add compilation-cache guidance and shape-change diagnostics so users can
   distinguish compilation time from solve time.
 
 ### 5.2 Make advanced geometry inference easier to set up safely
 
-- [ ] Accept named geometry/frame objects and prior specifications with units;
-  generate prior-predictive geometry plots and half-space checks before sampling.
+- [ ] Accept geometry keyword mappings, an explicit `LocalFrame`, and prior
+  mappings with units; generate prior-predictive geometry plots and half-space
+  checks before sampling.
 - [ ] Add multi-start geometry search and an initialization helper that can seed
   NUTS from deterministic fits while clearly separating optimization from
   posterior inference.
@@ -553,8 +570,9 @@ manually re-signing slip vectors.
   posterior output, and log-evidence estimates.
 - [ ] Validate SMC against analytic low-dimensional targets and NUTS on unimodal
   cases, then use it for deliberately multimodal geometry examples.
-- [ ] Define sampler-independent posterior/result protocols so BlackJAX API
-  changes or future samplers do not leak through the GeoDef user interface.
+- [ ] Define sampler-independent result records and sampling functions so
+  BlackJAX API changes or future samplers do not leak through the GeoDef user
+  interface.
 
 ---
 
@@ -567,11 +585,12 @@ does not multiply special cases in beginner-facing code.
 
 - [ ] Write a separate design note defining scope, state variables, sign/unit
   conventions, validation targets, and the boundary between static GeoDef
-  objects and a new optional `geodef.cycle` module.
+  values and a new optional `geodef.cycle` module.
 - [ ] Port and independently validate stress-kernel-driven quasi-dynamic
   rate-and-state evolution from `related/stress-shadows/unicycle/`.
-- [ ] Add friction-law objects, adaptive ODE integration, event detection,
-  restart/checkpoint files, and energy/moment diagnostics.
+- [ ] Add friction-law functions/callables, adaptive ODE integration, event
+  detection, restart/checkpoint files, and energy/moment diagnostics. Use a
+  state record only if checkpointing invariants require one.
 - [ ] Support rectangular and triangular faults, CPU first and differentiable
   JAX integration only after the reference CPU implementation is trusted.
 - [ ] Deliver a small pedagogical spring-slider example before a large fault-
@@ -584,7 +603,7 @@ does not multiply special cases in beginner-facing code.
 - [ ] Add layered half-space displacement Green's functions behind an optional
   dependency, beginning with a well-bounded elastic layering use case.
 - [ ] Evaluate viscoelastic and poroelastic engines only after source/engine
-  protocols can represent time and material parameters cleanly.
+  callable contracts can represent time and material parameters cleanly.
 - [ ] For every engine: cite equations, preserve a reference implementation,
   cross-validate published cases, declare capabilities/coordinate conventions,
   and show one end-to-end example through the same high-level workflow.
@@ -612,10 +631,12 @@ increments rather than becoming a long-lived rewrite.
 1. **v1.1.x consistency releases:** Priority 0, documentation corrections,
    packaging fixes, licensing/CI/typing scaffolding, cache-key completeness,
    validation, and deprecation scaffolding.
-2. **v1.2 beginner workflow:** `LocalFrame`, named geometry/slip/displacement,
-   dataset result views, friendly constructors, and the revised quickstart.
-3. **v1.3 problem and scale layer:** validated `SlipProblem`, noise operators,
-   nuisance parameters, module extractions, and large-problem diagnostics.
+2. **v1.2 beginner workflow:** a small function-oriented `invert`, `greens`,
+   `slip`, and `data` surface; `LocalFrame`; named result views; and the revised
+   quickstart.
+3. **v1.3 scale layer:** callable noise/linear operators, nuisance parameters,
+   module extractions, and large-problem diagnostics. `LinearSystem` remains the
+   optional prepared-system API for repeated analyses.
 4. **Parallel research releases:** remaining JAX work and SMC can proceed in
    small units once their touched public semantics are settled.
 5. **v2 candidates:** only genuinely breaking cleanup that survived a full
diff --git a/README.md b/README.md
index 2b04dfc..89b5bab 100644
--- a/README.md
+++ b/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 |
diff --git a/docs/backend.md b/docs/backend.md
index 02f2178..85a7238 100644
--- a/docs/backend.md
+++ b/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,
diff --git a/docs/bayes.md b/docs/bayes.md
index c3ea6b0..910fc5a 100644
--- a/docs/bayes.md
+++ b/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
 
diff --git a/docs/cache.md b/docs/cache.md
index a76031c..eddb496 100644
--- a/docs/cache.md
+++ b/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()
 ```
 
diff --git a/docs/conventions.md b/docs/conventions.md
index eb7a29e..b39369d 100644
--- a/docs/conventions.md
+++ b/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
 
diff --git a/docs/data.md b/docs/data.md
index 61ac997..1ad0915 100644
--- a/docs/data.md
+++ b/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
diff --git a/docs/fault.md b/docs/fault.md
index c2aa62b..93dc75f 100644
--- a/docs/fault.md
+++ b/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,9 +58,9 @@ 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`
@@ -62,7 +68,9 @@ Create a triangular fault directly from ENU vertex coordinates. Two forms:
   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
 ```
 
 ---
diff --git a/docs/geometry.md b/docs/geometry.md
new file mode 100644
index 0000000..dd47069
--- /dev/null
+++ b/docs/geometry.md
@@ -0,0 +1,103 @@
+# `geodef.geometry` — coordinate frames and array conversions
+
+> Conventions — axes, depth sign, angles, units, array ordering,
+> regularization: see [`conventions.md`](conventions.md).
+
+`Fault` is GeoDef's fault-geometry object. The `geometry` module supplies the
+smaller pieces needed to construct and transform its arrays without introducing
+a second geometry hierarchy.
+
+## `LocalFrame`
+
+`LocalFrame` gives local East/North/Up arrays a geographic origin. It is kept as
+a value object because the origin and projection must travel together whenever
+local arrays are stored, cached, or combined.
+
+```python
+import geodef
+
+frame = geodef.LocalFrame(
+    origin_lat=-2.0,
+    origin_lon=100.0,
+    origin_alt=0.0,
+)
+
+enu = frame.to_enu(lon=lon, lat=lat, alt=alt)
+geographic = frame.to_geographic(
+    east=enu[..., 0],
+    north=enu[..., 1],
+    up=enu[..., 2],
+)
+```
+
+`source.transform_enu(coordinates, target=target)` explicitly re-expresses an
+array in another frame. `source.require_compatible(other)` raises rather than
+silently combining coordinates defined by different origins.
+
+`Fault.planar(..., frame=frame)`, `Fault.from_triangles(..., frame=frame)`, and
+`Mesh(..., frame=frame)` attach the frame directly to the domain value that owns
+the coordinates.
+
+## Planar parameter vectors
+
+JAX geometry kernels use the compact expert vector
+`[e0, n0, depth, strike, dip, length, width]`. Deterministic and Bayesian
+geometry inference also accept a mapping, so callers need not rely on order:
+
+```python
+parameters = {
+    "e0": 0.0,
+    "n0": 0.0,
+    "depth": 15_000.0,
+    "strike": 315.0,
+    "dip": 25.0,
+    "length": 80_000.0,
+    "width": 40_000.0,
+}
+
+theta = geodef.geometry.as_planar_vector(parameters)
+parameters = geodef.geometry.planar_parameter_dict(theta)
+```
+
+Every value is validated for finiteness and physical range. The associated
+`LocalFrame` remains an explicit argument to geometry search or posterior
+construction.
+
+For ordinary forward models, construct the fault directly with named geographic
+keywords:
+
+```python
+fault = geodef.Fault.planar(
+    lat=-2.0,
+    lon=100.0,
+    depth=15_000.0,
+    strike=315.0,
+    dip=25.0,
+    length=80_000.0,
+    width=40_000.0,
+    n_length=12,
+    n_width=6,
+    frame=frame,
+)
+```
+
+## Triangular arrays
+
+Expand shared nodes and connectivity only when an engine needs per-triangle
+vertices:
+
+```python
+vertices = geodef.geometry.vertices_from_nodes(nodes_enu, triangles)
+strike, dip = geodef.geometry.triangle_strike_dip(vertices)
+
+fault = geodef.Fault.from_triangles(vertices, frame=frame)
+# Or avoid expansion at the call site:
+fault = geodef.Fault.from_triangles(
+    nodes_enu,
+    triangles=triangles,
+    frame=frame,
+)
+```
+
+`Fault.from_mesh(mesh)` is the shortest path when geographic nodes and
+connectivity already live in a `Mesh`.
diff --git a/docs/greens.md b/docs/greens.md
index 3e308d6..826a2a9 100644
--- a/docs/greens.md
+++ b/docs/greens.md
@@ -28,40 +28,47 @@ ordering as the stacked observation vector.
 
 ## High-level assembly
 
-### `greens(fault, datasets, *, components='both', rake=None, slip_azimuth=None) → np.ndarray`
+### `matrix(fault, datasets, *, components='both', rake=None, slip_azimuth=None, plate_rake=None) → np.ndarray`
 
 Build the projected Green's matrix for one or more datasets. Results are automatically cached.
 
 ```python
 import geodef
 
-G = geodef.greens.greens(fault, gnss)                # shape (n_obs, 2*N)
-G = geodef.greens.greens(fault, [gnss, insar])       # rows stacked vertically
+G = geodef.greens.matrix(fault, gnss)                # shape (n_obs, 2*N)
+G = geodef.greens.matrix(fault, [gnss, insar])       # rows stacked vertically
 ```
 
 By default columns are blocked: `[:N]` strike-slip, `[N:]` dip-slip.
 
-Pass `components=` to have `greens()` return a single-component matrix (shape
-`(n_obs, N)`) directly, using the same slip basis as `geodef.invert()`:
+Pass `components=` to have `matrix()` return a single-component matrix (shape
+`(n_obs, N)`) directly, using the same slip basis as `geodef.invert.solve()`:
 
 ```python
-G_strike = geodef.greens.greens(fault, gnss, components='strike')
-G_dip    = geodef.greens.greens(fault, gnss, components='dip')
-G_rake   = geodef.greens.greens(fault, gnss, components='rake', rake=90.0)
-G_az     = geodef.greens.greens(fault, gnss, components='azimuth', slip_azimuth=350.0)
+G_strike = geodef.greens.matrix(fault, gnss, components='strike')
+G_dip    = geodef.greens.matrix(fault, gnss, components='dip')
+G_rake   = geodef.greens.matrix(fault, gnss, components='rake', rake=90.0)
+G_az     = geodef.greens.matrix(fault, gnss, components='azimuth', slip_azimuth=350.0)
+G_plate  = geodef.greens.matrix(fault, gnss, components='plate', plate_rake=plate_rake)
 ```
 
 For `'rake'` (a single rake for every patch) and `'azimuth'` (a fixed
 geographic slip azimuth, so each patch's local rake is `slip_azimuth - strike_i`)
 the two blocked column sets are combined as `cos(theta)*G_strike + sin(theta)*G_dip`.
 
-### `select_slip_columns(G_full, n_patches, components, rake=None, fault_strike=None, slip_azimuth=None) → np.ndarray`
+`'plate'` keeps two blocked columns per patch, but rotates them to
+`[rake_parallel | rake_perpendicular]` using a scalar or per-patch
+`plate_rake`. This is the appropriate matrix basis when bounds and smoothing
+should follow a large-scale tectonic direction rather than variable local
+triangle orientations.
 
-The reduction primitive behind `greens(components=...)`. Apply it to any
+### `select_slip_columns(G_full, n_patches, components, rake=None, fault_strike=None, slip_azimuth=None, plate_rake=None) → np.ndarray`
+
+The reduction primitive behind `matrix(components=...)`. Apply it to any
 already-assembled `(M, 2*N)` matrix — a Green's matrix or a stress kernel — to
-project it into a one-component slip basis. `geodef.invert()` uses it internally
-for `components='strike'|'dip'|'rake'|'azimuth'` and to project stress-kernel
-regularization into the active basis.
+project it into a one-component slip basis. `geodef.invert.solve()` uses it internally
+for `components='strike'|'dip'|'rake'|'azimuth'|'plate'` and to project
+stress-kernel regularization into the active basis.
 
 ### `stack_obs(datasets) → np.ndarray`
 
@@ -79,7 +86,7 @@ Build block-diagonal inverse-covariance weight matrix.
 W = geodef.stack_weights([gnss, insar])   # shape (total_n_obs, total_n_obs)
 ```
 
-`geodef.invert()` calls these internally; use them directly when assembling `G @ m` by hand.
+`geodef.invert.solve()` calls these internally; use them directly when assembling `G @ m` by hand.
 
 `W` is the inverse covariance, not a vector of standard deviations. In
 derivations it is often clearer to whiten the system with a matrix `R` such
@@ -89,6 +96,17 @@ that `R.T @ R = W`, then solve with `R @ G` and `R @ d`.
 
 ## Laplacian operators
 
+### `project(data, G_raw) → np.ndarray`
+
+Project a raw three-component displacement matrix into a dataset's observed
+components. This is useful when assembling a custom forward operator.
+
+### `laplacian(fault) → np.ndarray`
+
+Return the fault's patch-order-aware Laplacian matrix. This is the convenient
+entry point for ordinary use; the builders below remain available for custom
+grids.
+
 ### `build_laplacian_2d(nL, nW) → np.ndarray`
 
 2-D finite-difference Laplacian for a structured rectangular grid. Each row sums to zero. Requires `nL >= 3` and `nW >= 3`.
@@ -146,14 +164,19 @@ Rectangular patches (Okada85). Returns shape `(3*nobs, 2*npatch)`.
 
 Rectangular patches, strain output. Shape `(4*nobs, 2*npatch)`. Pass `obs_depth` for internal points (uses Okada92).
 
-### `tri_displacement_greens(lat, lon, lat0, lon0, depth, vertices, nu=0.25)`
+### `tri_displacement_greens(lat, lon, lat0, lon0, depth, vertices, nu=0.25, *, frame=None)`
 
 Triangular patches (Nikkhoo & Walter). Returns shape `(3*nobs, 2*npatch)`.
 
-### `tri_strain_greens(lat, lon, lat0, lon0, depth, vertices, nu=0.25, obs_depth=None)`
+### `tri_strain_greens(lat, lon, lat0, lon0, depth, vertices, nu=0.25, obs_depth=None, *, frame=None)`
 
 Triangular patches, strain output. Shape `(6*nobs, 2*npatch)`.
 
+For direct low-level triangular calls, pass the `LocalFrame` that defines
+`vertices`. Omitting it retains legacy mean-centroid frame inference.
+`Fault.greens_matrix` always supplies `fault.frame`, so domain-level assembly
+cannot silently reinterpret triangular vertices.
+
 ---
 
 ## Column layout
diff --git a/docs/invert.md b/docs/invert.md
index 096e8be..28db637 100644
--- a/docs/invert.md
+++ b/docs/invert.md
@@ -28,24 +28,31 @@ for general background.
 
 ---
 
-## `invert(fault, datasets, **kwargs) → InversionResult`
+## `solve(fault, datasets, **kwargs) → InversionResult`
 
 ```python
 import geodef
 
 # Unregularized WLS
-result = geodef.invert(fault, [gnss, insar])
+result = geodef.invert.solve(fault, [gnss, insar])
 
 # Laplacian smoothing, non-negative
-result = geodef.invert(fault, [gnss, insar],
-                       smoothing='laplacian',
-                       smoothing_strength=1e3,
-                       bounds=(0, None))
+result = geodef.invert.solve(fault, [gnss, insar],
+                             smoothing='laplacian',
+                             smoothing_strength=1e3,
+                             bounds=(0, None))
 
 # One-parameter slip bases
-result = geodef.invert(fault, gnss, components='rake', rake=90.0)
-result = geodef.invert(fault, gnss,
-                       components='azimuth', slip_azimuth=15.0)
+result = geodef.invert.solve(fault, gnss, components='rake', rake=90.0)
+result = geodef.invert.solve(fault, gnss,
+                             components='azimuth', slip_azimuth=15.0)
+
+# Two plate-motion coordinates, suitable for variable-orientation meshes
+plate_rake = geodef.slip.plate_rake_from_euler(
+    fault, (pole_lat, pole_lon, rate)
+)
+result = geodef.invert.solve(fault, gnss,
+                             components='plate', plate_rake=plate_rake)
 ```
 
 ### Key parameters
@@ -55,11 +62,12 @@ result = geodef.invert(fault, gnss,
 | `method` | auto | `'wls'`, `'nnls'`, `'bounded_ls'`, `'constrained'` |
 | `smoothing` | `None` | `'laplacian'`, `'damping'`, `'stresskernel'`, or a custom matrix |
 | `smoothing_strength` | `0.0` | Regularization weight λ, or `'abic'`/`'cv'` for auto-tuning |
-| `smoothing_target` | `None` | Reference model for `(m - m_ref)` regularization |
+| `smoothing_target` | `None` | Vector reference for `(m - m_ref)` regularization |
 | `bounds` | `None` | `(lower, upper)` slip bounds; each side is a scalar, a per-component array, a per-parameter array, or `None` |
-| `components` | `'both'` | Slip basis: `'both'`, `'strike'`, `'dip'`, `'rake'`, or `'azimuth'` |
+| `components` | `'both'` | Slip basis: `'both'`, `'strike'`, `'dip'`, `'rake'`, `'azimuth'`, or `'plate'` |
 | `rake` | `None` | Fixed local rake angle in degrees; required for `components='rake'` |
 | `slip_azimuth` | `None` | Fixed geographic slip azimuth in degrees clockwise from north; required for `components='azimuth'` |
+| `plate_rake` | `None` | Scalar or per-patch large-scale direction in local rake coordinates; required for `components='plate'` |
 | `cv_folds` | `5` | Number of folds for cross-validation |
 | `constraints` | `None` | `(C, d)` for `C @ m <= d` (constrained solver only) |
 
@@ -86,14 +94,27 @@ using each patch's strike, so it is better for curved or variable-strike meshes.
 The Green's matrix and stress-kernel regularization are projected into the
 chosen slip basis automatically.
 
+`components='plate'` retains two parameters per patch but rotates them into
+large-scale rake-parallel/rake-perpendicular coordinates. Laplacian smoothing,
+targets, component bounds, covariance, and resolution all operate in that
+basis. This prevents abrupt triangle-local strike/dip changes from defining
+the regularization coordinates. The physical strike/dip slip remains available
+through `result.strike_slip` and `result.dip_slip`.
+
 ---
 
 ## `InversionResult`
 
 | Attribute | Shape | Description |
 |-----------|-------|-------------|
-| `slip` | `(N, 2)` or `(N, 1)` | Per-patch strike/dip slip for `components='both'`, or one active amplitude per patch |
-| `slip_vector` | `(2N,)` or `(N,)` | Blocked `[ss_0..ss_N, ds_0..ds_N]`, or one amplitude per patch |
+| `slip` | `(N, 2)` or `(N, 1)` | Backwards-compatible per-patch array in the solved coordinates |
+| `slip_vector` | `(2N,)` or `(N,)` | Backwards-compatible blocked vector in the solved coordinates |
+| `strike_slip` | `(N,)` | Physical strike-slip component |
+| `dip_slip` | `(N,)` | Physical dip-slip component |
+| `slip_magnitude` | `(N,)` | Unsigned physical slip magnitude |
+| `slip_rake` | `(N,)` | Physical local rake in degrees |
+| `rake_parallel` | `(N,)` | Plate-parallel solution block (`components='plate'` only) |
+| `rake_perpendicular` | `(N,)` | Plate-perpendicular solution block (`components='plate'` only) |
 | `predicted` | `(M,)` | Forward-modeled observations |
 | `residuals` | `(M,)` | `obs - predicted` |
 | `reduced_chi2` | scalar | Reduced chi-squared, `r^T W r / (M - P)` |
@@ -105,6 +126,10 @@ chosen slip basis automatically.
 | `components` | str | Slip basis used in the inversion |
 | `rake` | float or `None` | Fixed rake angle for `components='rake'` |
 | `slip_azimuth` | float or `None` | Fixed geographic azimuth for `components='azimuth'` |
+| `plate_rake` | `(N,)` or `None` | Per-patch large-scale direction for `components='plate'` |
+
+Use the named physical arrays for interpretation and plotting. Use
+`slip_vector` when assembling linear algebra in the solved basis.
 
 ```python
 result.save("result.npz")                   # save to disk
@@ -165,8 +190,8 @@ ac.optimal       # λ at minimum ABIC
 ### Auto-tuning via `smoothing_strength`
 
 ```python
-result = geodef.invert(fault, data, smoothing='laplacian', smoothing_strength='abic')
-result = geodef.invert(fault, data, smoothing='laplacian', smoothing_strength='cv')
+result = geodef.invert.solve(fault, data, smoothing='laplacian', smoothing_strength='abic')
+result = geodef.invert.solve(fault, data, smoothing='laplacian', smoothing_strength='cv')
 ```
 
 On the JAX backend (`geodef.backend.set_backend('jax')`), `abic_curve`
@@ -177,7 +202,7 @@ results, one fused sweep instead of a Python loop.
 
 ## Nonlinear geometry search (JAX)
 
-### `geometry_search(theta0, datasets, *, ref_lat, ref_lon, ...) → GeometrySearchResult`
+### `geometry_search(theta0, datasets, *, ...) → GeometrySearchResult`
 
 Gradient-based inversion for planar fault geometry: the slip is solved
 linearly inside a nonlinear search over selected geometry parameters
@@ -203,12 +228,15 @@ and `lambda` describe the chosen slip regularization.
 ```python
 geodef.backend.set_backend('jax')
 
-theta0 = [0.0, 0.0, 25e3, 315.0, 30.0, 180e3, 90e3]
-#         e0   n0   depth strike dip   length width  (start; true dip 15)
+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,
+}
 
 result = geodef.geometry_search(
-    theta0, gnss,
-    ref_lat=-2.0, ref_lon=100.0,     # anchors the local frame
+    geometry0, gnss, frame=frame,
     free=['dip', 'depth'],           # parameters to optimize; rest fixed
     bounds={'dip': (5.0, 45.0)},
     n_length=12, n_width=6,
@@ -216,7 +244,9 @@ result = geodef.geometry_search(
     smoothing='laplacian', smoothing_strength=1.0,
 )
 
-result.theta          # full 7-vector at the optimum
+result.fault          # concrete optimal Fault
+result.frame          # frame defining the local parameter vector
+result.theta          # expert/JAX seven-vector for the same geometry
 result.slip           # inner-solve slip at the optimal geometry
 result.theta_cov      # Gauss-Newton covariance of the free parameters
 result.reduced_chi2
@@ -224,8 +254,10 @@ result.reduced_chi2
 
 Notes:
 
-- `theta0` is in the local Cartesian frame anchored at
-  `(ref_lat, ref_lon)`; `e0`/`n0` are centroid offsets in meters.
+- For expert/JAX workflows, the seven-element `theta0` array remains
+  supported with either `frame=frame` or `ref_lat=..., ref_lon=...`.
+- `result.fault` is the ordinary domain view. `result.theta` is the exact
+  `[e0, n0, depth, strike, dip, length, width]` array view.
 - The inner solve is unconstrained WLS with fixed `smoothing_strength`;
   choose λ first (e.g. with `abic_curve` at a reasonable starting
   geometry).
diff --git a/docs/mesh.md b/docs/mesh.md
index 58cd556..9d256e0 100644
--- a/docs/mesh.md
+++ b/docs/mesh.md
@@ -29,11 +29,20 @@ from geodef.mesh import Mesh
 
 mesh.n_nodes       # number of vertices
 mesh.n_triangles   # number of triangles
+mesh.frame         # explicit LocalFrame, inferred from mean nodes by default
 mesh.centers_geo   # (M, 3) centroids as [lon, lat, depth_m]
 mesh.areas         # (M,) triangle areas in m²
-mesh.vertices_enu(ref_lat, ref_lon)  # (M, 3, 3) vertices in local ENU meters
+mesh.vertices_enu()                  # vertices in mesh.frame
+mesh.vertices_enu(frame=other_frame) # explicit alternate representation
 ```
 
+Pass `frame=geodef.LocalFrame(...)` when constructing a `Mesh` to choose the
+stored local representation. Legacy `vertices_enu(ref_lat, ref_lon)` remains
+supported. Supplying both forms is rejected because it would make provenance
+ambiguous.
+`mesh.to_frame(other_frame)` returns the same geographic mesh with a different
+default local representation.
+
 ### I/O
 
 ```python
diff --git a/docs/okada.md b/docs/okada.md
index 2216b72..66640e9 100644
--- a/docs/okada.md
+++ b/docs/okada.md
@@ -136,4 +136,4 @@ All low-level engines use a local Cartesian frame:
 - `z` — observation depth (≤ 0 for surface or below)
 - `depth` — fault centroid depth (positive down)
 
-The `Fault.greens_matrix()` and `greens.greens()` functions handle the geographic-to-local conversion automatically.
+The `Fault.greens_matrix()` and `greens.matrix()` functions handle the geographic-to-local conversion automatically.
diff --git a/docs/plot.md b/docs/plot.md
index 469a437..e9e89f3 100644
--- a/docs/plot.md
+++ b/docs/plot.md
@@ -25,7 +25,7 @@ interpolated map.
 Slip distribution as colored patches (rectangular or triangular).
 
 ```python
-geodef.plot.slip(fault, result.slip_vector)
+geodef.plot.slip(fault, result.slip_magnitude)
 
 geodef.plot.slip(fault, result.slip_vector,
     ax=ax,
@@ -52,9 +52,10 @@ station locations (the up-dip edge is then drawn as the surface trace). The
 `plot.resolution`, and `plot.uncertainty` (defaulting to `'geographic'` and
 `False` there).
 
-For one-parameter inversion results such as `components='rake'` or
-`components='azimuth'`, `result.slip_vector` has length `N`; `plot.slip()`
-plots that amplitude directly and ignores the `components` selector.
+Pass the named result array you want to plot, such as `result.strike_slip`,
+`result.dip_slip`, `result.slip_magnitude`, `result.rake_parallel`, or
+`result.rake_perpendicular`. Raw N/2N vectors remain accepted; a raw
+one-component vector is plotted directly.
 
 ---
 
@@ -66,8 +67,7 @@ Gouraud-shaded `pcolormesh` over the structured grid; triangular (or
 unstructured) faults use `tricontourf` over the patch centroids.
 
 ```python
-geodef.plot.slip_interpolated(fault, result.slip_vector,
-    components='magnitude',
+geodef.plot.slip_interpolated(fault, result.slip_magnitude,
     cmap='viridis',
     levels=20,          # filled contour levels (tricontourf path)
     colorbar=True,
diff --git a/docs/slip.md b/docs/slip.md
new file mode 100644
index 0000000..2c264df
--- /dev/null
+++ b/docs/slip.md
@@ -0,0 +1,107 @@
+# Slip vectors and basis conversions
+
+> Conventions — axes, depth sign, angles, units, array ordering,
+> regularization: see [`conventions.md`](conventions.md).
+
+GeoDef represents slip with NumPy arrays and provides functions for the common
+conversions. For `N` patches, a two-component model vector is blocked:
+`[strike_slip_0, ..., strike_slip_N, dip_slip_0, ..., dip_slip_N]`.
+
+## Pack and unpack
+
+```python
+from geodef import slip
+
+model = slip.pack(strike_slip, dip_slip)
+strike_slip, dip_slip = slip.unpack(model)
+
+slip_magnitude = slip.magnitude(strike_slip, dip_slip)
+slip_rake = slip.rake(strike_slip, dip_slip)
+```
+
+`fault.displacement` accepts the two physical components directly and returns
+three arrays:
+
+```python
+east, north, up = fault.displacement(
+    obs_lat,
+    obs_lon,
+    slip_strike=strike_slip,
+    slip_dip=dip_slip,
+)
+```
+
+## Fixed rake and geographic azimuth
+
+A one-component inversion solves one signed amplitude per patch. Convert an
+amplitude to physical strike/dip components with the matching function:
+
+```python
+strike_slip, dip_slip = slip.from_rake(amplitude, rake_degrees=90.0)
+
+strike_slip, dip_slip = slip.from_azimuth(
+    amplitude,
+    azimuth_degrees=15.0,
+    fault_strike_degrees=fault.strike,
+)
+```
+
+`from_rake` uses each patch's local strike/dip axes. `from_azimuth` preserves a
+single geographic direction across a curved mesh by accounting for each
+patch's strike.
+
+## Plate-motion coordinates
+
+For a curved triangular mesh, smoothing physical strike/dip components can
+inherit abrupt changes in patch orientation. A plate basis instead uses a
+smooth large-scale direction:
+
+```python
+plate_rake = slip.plate_rake_from_euler(
+    fault,
+    pole=(pole_lat, pole_lon, rate_degrees_per_myr),
+)
+
+strike_slip, dip_slip = slip.from_plate(
+    parallel,
+    perpendicular,
+    plate_rake_degrees=plate_rake,
+)
+
+parallel, perpendicular = slip.to_plate(
+    strike_slip,
+    dip_slip,
+    plate_rake_degrees=plate_rake,
+)
+```
+
+Invert directly in those coordinates with `components="plate"` and
+`plate_rake=plate_rake`. The result keeps the solved vector blocked as
+`[parallel | perpendicular]` and exposes both the solved and physical views:
+
+```python
+result = geodef.invert.solve(
+    fault,
+    datasets,
+    components="plate",
+    plate_rake=plate_rake,
+)
+
+result.rake_parallel
+result.rake_perpendicular
+result.strike_slip
+result.dip_slip
+result.slip_magnitude
+result.slip_rake
+```
+
+## Patch ordering
+
+Structured faults use `(n_width, n_length)` grid shape, with along-strike index
+varying fastest. Use the fault helpers instead of manual reshape assumptions:
+
+```python
+grid = fault.reshape_patches(per_patch_values)
+values = fault.flatten_patches(grid)
+index = fault.patch_index(strike_idx, dip_idx)
+```
diff --git a/examples/gorkha_earthquake/model_gorkha.ipynb b/examples/gorkha_earthquake/model_gorkha.ipynb
index 053d40a..d98b05c 100644
--- a/examples/gorkha_earthquake/model_gorkha.ipynb
+++ b/examples/gorkha_earthquake/model_gorkha.ipynb
@@ -140,7 +140,7 @@
     "# if we just accept the defaults, running the inversion is very easy:\n",
     "# this takes about 40 sec on my macbook the first time to compute the greens functions. \n",
     "# After the greens functions are cached, unsmoothed runs take < 1 sec, smoothed runs take about ~2sec.\n",
-    "result = geodef.invert(fault,[gnss,insar])\n",
+    "result = geodef.invert.solve(fault,[gnss,insar])\n",
     "\n",
     "# plot slip magnitude\n",
     "geodef.plot.slip(fault, result.slip_vector, coords=\"geographic\")\n",
@@ -196,7 +196,7 @@
    ],
    "source": [
     "# now that we have our optimum smoothing strength, we can run the smoothed inversion:\n",
-    "result = geodef.invert(fault,[gnss,insar],smoothing='laplacian',smoothing_strength=1e7)\n",
+    "result = geodef.invert.solve(fault,[gnss,insar],smoothing='laplacian',smoothing_strength=1e7)\n",
     "\n",
     "# plot slip magnitude\n",
     "ax=geodef.plot.slip(fault, result.slip_vector, coords=\"geographic\")\n",
@@ -299,7 +299,7 @@
    ],
    "source": [
     "# Now fit a non-negative model with fixed geographic slip azimuth.\n",
-    "result = geodef.invert(\n",
+    "result = geodef.invert.solve(\n",
     "    fault, [gnss, insar],\n",
     "    smoothing='laplacian', smoothing_strength=1e8,\n",
     "    components='azimuth', slip_azimuth=15.0,\n",
diff --git a/src/geodef/__init__.py b/src/geodef/__init__.py
index 45d287c..232ce17 100644
--- a/src/geodef/__init__.py
+++ b/src/geodef/__init__.py
@@ -8,20 +8,24 @@
     cache,
     euler,
     geomap,
+    geometry,
     gradients,
     greens,
+    invert,
     medium,
     mesh,
     okada,
     okada85,
     okada92,
     plot,
+    slip,
     transforms,
     tri,
     validation,
 )
 from geodef.data import GNSS, DataSet, InSAR, Vertical, spatial_covariance
 from geodef.fault import Fault, magnitude_to_moment, moment_to_magnitude
+from geodef.geometry import LocalFrame
 from geodef.greens import select_slip_columns, stack_obs, stack_weights
 from geodef.invert import (
     ABICCurveResult,
@@ -34,11 +38,11 @@
     compute_abic,
     dataset_diagnostics,
     geometry_search,
-    invert,
     lcurve,
     model_covariance,
     model_resolution,
     model_uncertainty,
+    solve,
 )
 from geodef.medium import DEFAULT_MEDIUM, ElasticMedium
 
@@ -49,14 +53,17 @@
     "cache",
     "euler",
     "geomap",
+    "geometry",
     "gradients",
     "greens",
+    "invert",
     "medium",
     "mesh",
     "okada",
     "okada85",
     "okada92",
     "plot",
+    "slip",
     "transforms",
     "tri",
     "validation",
@@ -68,6 +75,7 @@
     "spatial_covariance",
     # Fault geometry, medium, and moment
     "Fault",
+    "LocalFrame",
     "ElasticMedium",
     "DEFAULT_MEDIUM",
     "magnitude_to_moment",
@@ -87,7 +95,7 @@
     "compute_abic",
     "dataset_diagnostics",
     "geometry_search",
-    "invert",
+    "solve",
     "lcurve",
     "model_covariance",
     "model_resolution",
diff --git a/src/geodef/bayes.py b/src/geodef/bayes.py
index f23aa9e..621f355 100644
--- a/src/geodef/bayes.py
+++ b/src/geodef/bayes.py
@@ -30,8 +30,9 @@
 
     import geodef
     geodef.backend.set_backend("jax")
+    frame = geodef.LocalFrame(-2.0, 100.0)
     post = geodef.bayes.RectPosterior(
-        theta0, datasets, ref_lat=..., ref_lon=...,
+        theta0, datasets, frame=frame,
         free=["dip", "depth"],
         theta_prior={"dip": (5.0, 60.0), "depth": (5e3, 40e3)},
         n_length=8, n_width=4, smoothing="laplacian",
@@ -42,20 +43,27 @@
 from __future__ import annotations
 
 import dataclasses
-from collections.abc import Sequence
+from collections.abc import Mapping, Sequence
 from typing import Any, cast
 
 import numpy as np
 import numpy.typing as npt
 import scipy.linalg
 
-from geodef import backend, transforms
+from geodef import backend
 from geodef.data import DataSet
 from geodef.fault import Fault
+from geodef.geometry import (
+    LocalFrame,
+    _resolve_frame,
+    as_planar_vector,
+    planar_parameter_dict,
+)
 from geodef.gradients import rect_greens, tri_greens
 from geodef.invert import (
     _THETA_NAMES,
     LinearSystem,
+    _fault_from_planar_vector,
     _projection_matrix,
     _rank_positive_eigs,
 )
@@ -552,14 +560,16 @@ class RectPosterior(_CollapsedPosterior):
     ``-inf``.
 
     Args:
-        theta0: Template geometry ``[e0, n0, depth, strike, dip,
-            length, width]``; fixed parameters keep these values and
-            free ones are initialized from them. ``e0``/``n0`` are
-            centroid offsets in meters from (ref_lat, ref_lon).
+        theta0: Named parameter mapping, or expert array
+            ``[e0, n0, depth, strike, dip, length, width]``. Fixed parameters
+            keep these values and free ones are initialized from them. Array
+            input requires ``frame`` or ``ref_lat``/``ref_lon``.
         datasets: One or more displacement datasets (GNSS, InSAR,
             Vertical).
         ref_lat: Latitude anchoring the local Cartesian frame.
         ref_lon: Longitude anchoring the local Cartesian frame.
+        frame: Explicit local frame for array ``theta0``. Mutually exclusive
+            with an incompatible legacy ``ref_lat``/``ref_lon`` origin.
         free: Names of geometry parameters to sample. May be empty for
             pure hyperparameter inference.
         theta_prior: Prior for each free geometry parameter, keyed by
@@ -607,11 +617,12 @@ class RectPosterior(_CollapsedPosterior):
 
     def __init__(
         self,
-        theta0: npt.ArrayLike,
+        theta0: npt.ArrayLike | Mapping[str, float],
         datasets: DataSet | list[DataSet],
         *,
-        ref_lat: float,
-        ref_lon: float,
+        ref_lat: float | None = None,
+        ref_lon: float | None = None,
+        frame: LocalFrame | None = None,
         free: Sequence[str] = ("depth", "dip"),
         theta_prior: dict[str, tuple] | None = None,
         n_length: int = 1,
@@ -659,12 +670,15 @@ def __init__(
                 "mode='profiled' requires a fixed smoothing_strength (lambda)"
             )
 
-        theta0 = np.asarray(theta0, dtype=float)
+        frame = _resolve_frame(frame, ref_lat, ref_lon)
+        theta0 = as_planar_vector(theta0)
         self.mode = mode
         self.free = list(free)
         self.components = components
         self._components = components
         self.datasets = datasets
+        self.theta0 = np.array(theta0, copy=True)
+        self.frame = frame
         self._theta0 = theta0
         self._free_idx = np.array(
             [_THETA_NAMES.index(name) for name in free], dtype=int
@@ -675,17 +689,7 @@ def __init__(
 
         # Template system provides the stacked data, weights, and
         # regularization operator; its Green's matrix is not used.
-        template = Fault.planar(
-            lat=ref_lat,
-            lon=ref_lon,
-            depth=theta0[2],
-            strike=theta0[3],
-            dip=theta0[4],
-            length=theta0[5],
-            width=theta0[6],
-            n_length=n_length,
-            n_width=n_width,
-        )
+        template = _fault_from_planar_vector(theta0, frame, n_length, n_width)
         sys = LinearSystem(template, datasets, smoothing, components)
         n_patches = n_length * n_width
         self._col_start, self._col_stop = {
@@ -699,11 +703,13 @@ def __init__(
 
         e_parts, n_parts = [], []
         for ds in datasets:
-            e_ds, n_ds, _ = transforms.geod2enu(
-                ds.lat, ds.lon, np.zeros(ds.n_stations), ref_lat, ref_lon, 0.0
+            enu = frame.to_enu(
+                lon=ds.lon,
+                lat=ds.lat,
+                alt=np.full(ds.n_stations, frame.origin_alt),
             )
-            e_parts.append(e_ds)
-            n_parts.append(n_ds)
+            e_parts.append(enu[:, 0])
+            n_parts.append(enu[:, 1])
         self._e_obs = np.concatenate(e_parts)
         self._n_obs = np.concatenate(n_parts)
 
@@ -787,6 +793,42 @@ def __init__(
             self._mask = np.zeros(n_params, dtype=bool)
             self._logpdf_fn = self._build_logpdf()
 
+    def geometry(self, x: npt.ArrayLike) -> dict[str, float]:
+        """Return named planar parameters represented by one state.
+
+        This is a user-facing, non-JAX view. The likelihood methods continue to
+        consume array states directly for tracing and vectorization.
+
+        Args:
+            x: One posterior parameter state, shape ``(n_params,)``.
+
+        Returns:
+            Parameter dictionary in the posterior's :attr:`frame`.
+
+        Raises:
+            ValueError: If ``x`` is not one complete parameter state.
+        """
+        state = np.asarray(x, dtype=float)
+        if state.shape != (self.n_params,):
+            raise ValueError(f"x must have shape ({self.n_params},), got {state.shape}")
+        theta = np.array(self._theta0, copy=True)
+        theta[self._free_idx] = state[: len(self.free)]
+        return planar_parameter_dict(theta)
+
+    def fault(self, x: npt.ArrayLike) -> Fault:
+        """Return the planar fault represented by one parameter state.
+
+        Args:
+            x: One posterior parameter state, shape ``(n_params,)``.
+
+        Returns:
+            A fault discretized with this posterior's grid shape.
+        """
+        theta = as_planar_vector(self.geometry(x))
+        return _fault_from_planar_vector(
+            theta, self.frame, self._n_length, self._n_width
+        )
+
     def _setup_positive(
         self,
         positive: str | npt.ArrayLike,
@@ -1104,8 +1146,8 @@ class TriWarp:
     the same interpolated offset, since they share the same (u, v)).
 
     Args:
-        fault: Reference triangular ``Fault`` (``fault.vertices`` must not
-            be None); kept for :meth:`fault` and by :class:`TriPosterior`.
+        fault: Reference triangular :class:`Fault`; kept for :meth:`fault` and
+            by :class:`TriPosterior`.
         knots: Explicit knot locations in the mesh's best-fit-plane (u, v)
             coordinates, shape (nk, 2). Takes precedence over ``n_knots``.
         n_knots: ``(n_u, n_v)`` grid shape spanning the mesh's (u, v)
@@ -1118,9 +1160,8 @@ class TriWarp:
             solving for the interpolation weights (numerical stability).
 
     Raises:
-        ValueError: If ``fault`` is not a triangular fault, ``knots`` has
-            the wrong shape, or no default ``length_scale`` can be
-            inferred.
+        ValueError: If ``fault`` is not triangular, ``knots`` has the wrong
+            shape, or no default ``length_scale`` can be inferred.
     """
 
     def __init__(
@@ -1137,8 +1178,9 @@ def __init__(
                 "TriWarp requires a triangular Fault (fault.vertices is None)"
             )
         self._ref_fault = fault
-        self._ref_lat = float(np.mean(fault.centers[:, 0]))
-        self._ref_lon = float(np.mean(fault.centers[:, 1]))
+        self.frame = fault.frame
+        self._ref_lat = self.frame.origin_lat
+        self._ref_lon = self.frame.origin_lon
 
         v0 = np.asarray(fault.vertices, dtype=float)
         self._shape = v0.shape
@@ -1279,7 +1321,9 @@ def fault(self, theta: npt.ArrayLike) -> Fault:
         """
         verts = backend.to_numpy(self.vertices(np.asarray(theta, dtype=float)))
         return Fault.from_triangles(
-            verts.astype(float), ref_lat=self._ref_lat, ref_lon=self._ref_lon
+            verts.astype(float),
+            frame=self.frame,
+            medium=self._ref_fault.medium,
         )
 
     def plot(self, theta: npt.ArrayLike | None = None, ax: Any = None) -> tuple:
@@ -1494,14 +1538,17 @@ def __init__(
         self._n_slip = n_slip
         self._n_patches = n_patches
 
-        ref_lat, ref_lon = warp._ref_lat, warp._ref_lon
+        frame = warp.frame
+        self.frame = frame
         e_parts, n_parts = [], []
         for ds in datasets:
-            e_ds, n_ds, _ = transforms.geod2enu(
-                ds.lat, ds.lon, np.zeros(ds.n_stations), ref_lat, ref_lon, 0.0
+            enu = frame.to_enu(
+                lon=ds.lon,
+                lat=ds.lat,
+                alt=np.full(ds.n_stations, frame.origin_alt),
             )
-            e_parts.append(e_ds)
-            n_parts.append(n_ds)
+            e_parts.append(enu[:, 0])
+            n_parts.append(enu[:, 1])
         e_obs = np.concatenate(e_parts)
         n_obs = np.concatenate(n_parts)
         self._obs = np.column_stack([e_obs, n_obs, np.zeros_like(e_obs)])
diff --git a/src/geodef/data.py b/src/geodef/data.py
index f80ea72..5235e89 100644
--- a/src/geodef/data.py
+++ b/src/geodef/data.py
@@ -861,7 +861,7 @@ def spatial_covariance(
     This is the standard way to represent spatially-correlated InSAR noise
     (atmosphere, orbits) instead of assuming diagonal ``C_d``. Pass the result
     as ``covariance=`` to a single-component-per-station dataset (``InSAR`` or
-    ``Vertical``), or thread it through ``geodef.invert()``.
+    ``Vertical``), or thread it through ``geodef.invert.solve()``.
 
     The covariance is ``C_ij = sill * rho(d_ij) + nugget * delta_ij`` where
     ``rho`` is the correlation function:
diff --git a/src/geodef/fault.py b/src/geodef/fault.py
index 2d2e2ba..b19191d 100644
--- a/src/geodef/fault.py
+++ b/src/geodef/fault.py
@@ -8,9 +8,16 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+import numpy.typing as npt
 
 from geodef import greens as _greens
 from geodef import transforms
+from geodef.geometry import (
+    LocalFrame,
+    as_planar_vector,
+    triangle_strike_dip,
+    vertices_from_nodes,
+)
 from geodef.medium import DEFAULT_MEDIUM, ElasticMedium
 from geodef.validation import (
     ValidationReport,
@@ -48,6 +55,8 @@ class Fault:
         medium: Elastic half-space parameters used by Green's functions,
             stress kernels, and moment. Defaults to
             ``geodef.medium.DEFAULT_MEDIUM`` (30 GPa Poisson solid).
+        frame: Local frame for local-coordinate views. Inferred from mean patch
+            coordinates for direct legacy construction.
     """
 
     def __init__(
@@ -64,6 +73,7 @@ def __init__(
         grid_shape: tuple[int, int] | None = None,
         engine: str = "okada",
         medium: ElasticMedium | None = None,
+        frame: LocalFrame | None = None,
     ) -> None:
         lat = as_1d_floats("lat", np.atleast_1d(lat), unit="degrees")
         n = lat.shape[0]
@@ -122,6 +132,8 @@ def __init__(
         self._grid_shape = grid_shape
         self._engine = engine
         self._medium = DEFAULT_MEDIUM if medium is None else medium
+        inferred_frame = LocalFrame(float(np.mean(lat)), float(np.mean(lon)))
+        self._frame = inferred_frame if frame is None else frame
 
         # Make arrays read-only
         for arr in (self._lat, self._lon, self._depth, self.strike, self.dip):
@@ -135,8 +147,8 @@ def __init__(
         # Lazy caches
         self._centers_local: np.ndarray | None = None
         self._laplacian: np.ndarray | None = None
-        self._ref_lat = float(np.mean(lat))
-        self._ref_lon = float(np.mean(lon))
+        self._ref_lat = self._frame.origin_lat
+        self._ref_lon = self._frame.origin_lon
 
     # ==================================================================
     # Factory classmethods
@@ -156,6 +168,7 @@ def planar(
         n_length: int = 1,
         n_width: int = 1,
         medium: ElasticMedium | None = None,
+        frame: LocalFrame | None = None,
     ) -> "Fault":
         """Create a discretized planar fault from its center.
 
@@ -171,20 +184,17 @@ def planar(
             n_width: Number of patches down dip.
             medium: Elastic half-space parameters. Defaults to the 30 GPa
                 Poisson solid ``geodef.medium.DEFAULT_MEDIUM``.
+            frame: Local frame for local-coordinate views. Defaults to a frame
+                centered horizontally on the fault.
 
         Returns:
             A Fault with ``n_length * n_width`` rectangular patches.
         """
-        for pname, val, unit in (
-            ("lat", lat, "degrees"),
-            ("lon", lon, "degrees"),
-            ("depth", depth, "meters"),
-            ("strike", strike, "degrees"),
-            ("dip", dip, "degrees"),
-            ("length", length, "meters"),
-            ("width", width, "meters"),
-        ):
-            check_finite_scalar(pname, val, unit=unit)
+        if n_length < 1 or n_width < 1:
+            raise ValueError("n_length and n_width must be positive integers")
+        selected_frame = LocalFrame(lat, lon) if frame is None else frame
+        center = selected_frame.to_enu(lon=lon, lat=lat, alt=selected_frame.origin_alt)
+        as_planar_vector([center[0], center[1], depth, strike, dip, length, width])
         patch_L = length / n_length
         patch_W = width / n_width
 
@@ -217,14 +227,13 @@ def planar(
         )
         u_offsets = fault_u0 + (jj + 0.5) * patch_W * sin_dip
 
-        lat_c, lon_c, _ = transforms.translate_flat(
-            lat,
-            lon,
-            0.0,
-            e_offsets,
-            n_offsets,
-            0.0,
+        geographic = selected_frame.to_geographic(
+            east=center[0] + e_offsets,
+            north=center[1] + n_offsets,
+            up=np.zeros_like(e_offsets),
         )
+        lon_c = geographic[:, 0]
+        lat_c = geographic[:, 1]
         depth_c = depth - u_offsets
 
         n_patches = n_length * n_width
@@ -244,6 +253,7 @@ def planar(
             grid_shape=(n_length, n_width),
             engine="okada",
             medium=medium,
+            frame=selected_frame,
         )
 
     @classmethod
@@ -330,8 +340,9 @@ def from_triangles(
         cls,
         vertices: np.ndarray,
         *,
-        ref_lat: float = 0.0,
-        ref_lon: float = 0.0,
+        ref_lat: float | None = None,
+        ref_lon: float | None = None,
+        frame: LocalFrame | None = None,
         triangles: np.ndarray | None = None,
         medium: ElasticMedium | None = None,
     ) -> "Fault":
@@ -350,10 +361,13 @@ def from_triangles(
           mesh and its node sharing.
 
         Args:
-            vertices: Either per-triangle corners, shape (N, 3, 3), or a shared
-                node array, shape (M, 3), each row [east, north, up] in meters.
+            vertices: Per-triangle corners with shape (N, 3, 3), or a shared
+                node array with shape (M, 3). Numeric rows are [east, north,
+                up] in meters.
             ref_lat: Reference latitude for the ENU origin.
             ref_lon: Reference longitude for the ENU origin.
+            frame: Explicit local frame. Mutually exclusive with legacy
+                ``ref_lat``/``ref_lon``.
             triangles: Optional connectivity indices into ``vertices``, shape
                 (N, 3). When given, ``vertices`` is treated as a node array.
             medium: Elastic half-space parameters. Defaults to the 30 GPa
@@ -366,50 +380,48 @@ def from_triangles(
             ValueError: If the array shapes are inconsistent, or a triangle
                 index is out of range.
         """
-        from geodef.mesh import _compute_strike_dip
-
-        vertices = np.asarray(vertices, dtype=float)
-        if triangles is not None:
-            triangles = np.asarray(triangles, dtype=int)
-            if vertices.ndim != 2 or vertices.shape[1] != 3:
+        legacy_frame: LocalFrame | None = None
+        if (ref_lat is None) != (ref_lon is None):
+            raise ValueError("ref_lat and ref_lon must be provided together")
+        if ref_lat is not None and ref_lon is not None:
+            legacy_frame = LocalFrame(ref_lat, ref_lon)
+        if frame is not None and legacy_frame is not None:
+            frame.require_compatible(legacy_frame)
+
+        selected_frame = frame if frame is not None else legacy_frame
+        if selected_frame is None:
+            selected_frame = LocalFrame(0.0, 0.0)
+        if triangles is None:
+            vertices_array = np.asarray(vertices, dtype=float)
+            if vertices_array.ndim != 3 or vertices_array.shape[1:] != (3, 3):
                 raise ValueError(
-                    "with triangles, vertices must be a node array of shape (M, 3)"
+                    "without triangles, vertices must have shape (N, 3, 3)"
                 )
-            if triangles.ndim != 2 or triangles.shape[1] != 3:
-                raise ValueError("triangles must have shape (N, 3)")
-            if triangles.size and (
-                triangles.min() < 0 or triangles.max() >= vertices.shape[0]
-            ):
-                raise ValueError("triangles index out of range for the node array")
-            vertices = vertices[triangles]  # (N, 3, 3)
-        elif vertices.ndim != 3 or vertices.shape[1:] != (3, 3):
-            raise ValueError("vertices must have shape (N, 3, 3)")
-
-        strike, dip = _compute_strike_dip(vertices)
-
-        # Compute centroids in ENU then convert to geographic
-        centroids_enu = np.mean(vertices, axis=1)  # (N, 3)
-        lat, lon, alt = transforms.enu2geod(
-            centroids_enu[:, 0],
-            centroids_enu[:, 1],
-            centroids_enu[:, 2],
-            ref_lat,
-            ref_lon,
-            0.0,
+            if not np.all(np.isfinite(vertices_array)):
+                raise ValueError("vertices must contain only finite values")
+        else:
+            vertices_array = vertices_from_nodes(vertices, triangles)
+
+        centers_enu = np.mean(vertices_array, axis=1)
+        centers_geo = selected_frame.to_geographic(
+            east=centers_enu[:, 0],
+            north=centers_enu[:, 1],
+            up=centers_enu[:, 2],
         )
-        depth = -alt  # ENU up → depth positive down
+        strike, dip = triangle_strike_dip(vertices_array)
 
         return cls(
-            lat,
-            lon,
-            depth,
+            centers_geo[:, 1],
+            centers_geo[:, 0],
+            -centers_geo[:, 2],
             strike,
             dip,
             None,
             None,
-            vertices=vertices,
+            vertices=vertices_array,
             engine="tri",
             medium=medium,
+            frame=selected_frame,
         )
 
     @classmethod
@@ -426,15 +438,25 @@ def from_mesh(
 
         Args:
             mesh: A ``geodef.mesh.Mesh`` instance.
+            medium: Elastic half-space parameters. Defaults to the 30 GPa
+                Poisson solid ``geodef.medium.DEFAULT_MEDIUM``.
 
         Returns:
             A triangular Fault with ``engine="tri"``.
         """
-        ref_lat = float(np.mean(mesh.lat))
-        ref_lon = float(np.mean(mesh.lon))
-        vertices = mesh.vertices_enu(ref_lat, ref_lon)
+        frame = mesh.frame
+        if frame is None:
+            raise ValueError("mesh must define a local frame")
+        nodes_enu = frame.to_enu(
+            lon=mesh.lon,
+            lat=mesh.lat,
+            alt=-mesh.depth,
+        )
         return cls.from_triangles(
-            vertices, ref_lat=ref_lat, ref_lon=ref_lon, medium=medium
+            nodes_enu,
+            triangles=mesh.triangles,
+            frame=frame,
+            medium=medium,
         )
 
     @classmethod
@@ -733,21 +755,27 @@ def centers_geo(self) -> np.ndarray:
     def centers_local(self) -> np.ndarray:
         """Patch centroids in local Cartesian [east, north, up] in meters.
 
-        Computed relative to the fault centroid (mean lat/lon).
+        Coordinates are expressed in :attr:`frame`.
         """
         if self._centers_local is None:
-            alt = np.zeros(self.n_patches)
-            e, n, u = transforms.geod2enu(
-                self._lat,
-                self._lon,
-                alt,
-                self._ref_lat,
-                self._ref_lon,
-                0.0,
-            )
-            self._centers_local = np.column_stack([e, n, -self._depth])
+            if self._vertices is not None:
+                self._centers_local = np.mean(self._vertices, axis=1)
+            else:
+                enu = self._frame.to_enu(
+                    lon=self._lon,
+                    lat=self._lat,
+                    alt=np.full(self.n_patches, self._frame.origin_alt),
+                )
+                self._centers_local = np.column_stack(
+                    [enu[:, 0], enu[:, 1], -self._depth]
+                )
         return self._centers_local
 
+    @property
+    def frame(self) -> LocalFrame:
+        """Local coordinate frame for :attr:`centers_local` and vertices."""
+        return self._frame
+
     @property
     def areas(self) -> np.ndarray:
         """Patch areas in square meters, shape (N,)."""
@@ -852,6 +880,40 @@ def with_medium(self, medium: ElasticMedium) -> "Fault":
             grid_shape=self._grid_shape,
             engine=self._engine,
             medium=medium,
+            frame=self._frame,
+        )
+
+    def to_frame(self, frame: LocalFrame) -> "Fault":
+        """Return this fault explicitly re-expressed in another local frame.
+
+        Geographic rectangular patch coordinates remain unchanged. Triangular
+        vertices are transformed so their physical geographic positions remain
+        unchanged rather than being reinterpreted in the new frame.
+
+        Args:
+            frame: Destination local frame.
+
+        Returns:
+            A fault with the same physical geometry and elastic medium in
+            ``frame``.
+        """
+        if self._frame.is_compatible(frame):
+            return self
+        if self._vertices is not None:
+            vertices = self._frame.transform_enu(self._vertices, target=frame)
+            return Fault.from_triangles(vertices, frame=frame, medium=self._medium)
+        return Fault(
+            self._lat,
+            self._lon,
+            self._depth,
+            self.strike,
+            self.dip,
+            self._length,
+            self._width,
+            grid_shape=self._grid_shape,
+            engine=self._engine,
+            medium=self._medium,
+            frame=frame,
         )
 
     @property
@@ -961,6 +1023,7 @@ def greens_matrix(
                     self._depth,
                     self._vertices,
                     nu=nu,
+                    frame=self._frame,
                 )
             elif kind == "strain":
                 return _greens.tri_strain_greens(
@@ -972,6 +1035,7 @@ def greens_matrix(
                     self._vertices,
                     nu=nu,
                     obs_depth=obs_depth,
+                    frame=self._frame,
                 )
             raise ValueError(f"Unknown kind: {kind!r}. Use 'displacement' or 'strain'.")
 
@@ -989,12 +1053,12 @@ def displacement(
         Args:
             obs_lat: Observation latitudes, shape (M,).
             obs_lon: Observation longitudes, shape (M,).
-            slip_strike: Strike-slip component per patch. Scalar (broadcast
-                to all patches) or array of shape (N,).
+            slip_strike: Strike-slip component per patch. Values may be scalar
+                (broadcast to all patches) or shape (N,).
             slip_dip: Dip-slip component per patch. Scalar or array of shape (N,).
 
         Returns:
-            Tuple (ue, un, uz) of displacement arrays, each shape (M,).
+            ``(east, north, up)`` displacement arrays, each shape (M,).
         """
         obs_lat = np.atleast_1d(np.asarray(obs_lat, dtype=float))
         obs_lon = np.atleast_1d(np.asarray(obs_lon, dtype=float))
@@ -1057,7 +1121,7 @@ def stress_kernel(self, mu: float | None = None) -> np.ndarray:
     # Moment and magnitude
     # ==================================================================
 
-    def moment(self, slip: np.ndarray, mu: float | None = None) -> float:
+    def moment(self, slip: npt.ArrayLike, mu: float | None = None) -> float:
         """Compute scalar seismic moment.
 
         Args:
@@ -1070,10 +1134,10 @@ def moment(self, slip: np.ndarray, mu: float | None = None) -> float:
         """
         if mu is None:
             mu = self._medium.shear_modulus
-        slip = np.asarray(slip, dtype=float)
-        return float(mu * np.sum(slip * self.areas))
+        slip_array = np.asarray(slip, dtype=float)
+        return float(mu * np.sum(slip_array * self.areas))
 
-    def magnitude(self, slip: np.ndarray, mu: float | None = None) -> float:
+    def magnitude(self, slip: npt.ArrayLike, mu: float | None = None) -> float:
         """Compute moment magnitude from a slip distribution.
 
         Args:
@@ -1111,6 +1175,58 @@ def patch_index(self, strike_idx: int, dip_idx: int) -> int:
         nL, _ = self._grid_shape
         return dip_idx * nL + strike_idx
 
+    def reshape_patches(self, values: npt.ArrayLike) -> np.ndarray:
+        """Reshape a patch-first array into ``[dip, strike, ...]`` grid axes.
+
+        The helper makes the storage convention explicit: strike index varies
+        fastest, so a flat vector becomes a grid of shape
+        ``(n_width, n_length)``. Trailing value dimensions are preserved.
+
+        Args:
+            values: Array whose first axis has length ``n_patches``.
+
+        Returns:
+            Array with leading axes ``(n_width, n_length)``.
+
+        Raises:
+            ValueError: If the fault is unstructured or the leading dimension
+                does not match the patch count.
+        """
+        if self._grid_shape is None:
+            raise ValueError("reshape_patches requires a structured grid")
+        array = np.asarray(values)
+        if array.ndim == 0 or array.shape[0] != self.n_patches:
+            raise ValueError(
+                f"values must have leading dimension {self.n_patches}, got "
+                f"shape {array.shape}"
+            )
+        n_length, n_width = self._grid_shape
+        return array.reshape((n_width, n_length, *array.shape[1:]))
+
+    def flatten_patches(self, values: npt.ArrayLike) -> np.ndarray:
+        """Flatten ``[dip, strike, ...]`` grid axes into patch storage order.
+
+        Args:
+            values: Array with leading shape ``(n_width, n_length)``.
+
+        Returns:
+            Patch-first array with leading dimension ``n_patches``.
+
+        Raises:
+            ValueError: If the fault is unstructured or leading grid axes do
+                not match the fault.
+        """
+        if self._grid_shape is None:
+            raise ValueError("flatten_patches requires a structured grid")
+        array = np.asarray(values)
+        n_length, n_width = self._grid_shape
+        expected = (n_width, n_length)
+        if array.ndim < 2 or array.shape[:2] != expected:
+            raise ValueError(
+                f"values must have leading shape {expected}, got shape {array.shape}"
+            )
+        return array.reshape((self.n_patches, *array.shape[2:]))
+
     # ==================================================================
     # File I/O
     # ==================================================================
@@ -1315,16 +1431,14 @@ def _save_tri_ned(self, fname: str) -> None:
         n_tri = self.n_patches
         verts_flat = self._vertices.reshape(-1, 3)  # (N*3, 3) [east, north, up]
 
-        # Convert ENU offsets (relative to fault centroid) back to geographic
-        lat_nodes, lon_nodes, _ = transforms.translate_flat(
-            self._ref_lat,
-            self._ref_lon,
-            0.0,
-            verts_flat[:, 0],
-            verts_flat[:, 1],
-            0.0,
+        geographic = self._frame.to_geographic(
+            east=verts_flat[:, 0],
+            north=verts_flat[:, 1],
+            up=verts_flat[:, 2],
         )
-        depth_nodes = -verts_flat[:, 2]  # up -> positive-down depth
+        lon_nodes = geographic[:, 0]
+        lat_nodes = geographic[:, 1]
+        depth_nodes = -geographic[:, 2]
 
         # Deduplicate nodes with fixed precision to merge shared vertices
         coords = np.column_stack([lon_nodes, lat_nodes, depth_nodes])
@@ -1338,6 +1452,7 @@ def _save_tri_ned(self, fname: str) -> None:
             lat=lat_nodes[unique_idx],
             depth=depth_nodes[unique_idx],
             triangles=inverse.reshape(n_tri, 3),
+            frame=self._frame,
         )
         mesh.save(fname)
 
@@ -1376,16 +1491,17 @@ def to_gmt(
             verts_enu = self._vertices  # (N, 3, 3)
             assert verts_enu is not None
             verts_flat = verts_enu.reshape(-1, 3)
-            lat_v, lon_v, _ = transforms.translate_flat(
-                self._ref_lat,
-                self._ref_lon,
-                0.0,
-                verts_flat[:, 0],
-                verts_flat[:, 1],
-                0.0,
+            geographic = self._frame.to_geographic(
+                east=verts_flat[:, 0],
+                north=verts_flat[:, 1],
+                up=verts_flat[:, 2],
             )
             verts = np.stack(
-                [lon_v.reshape(n, 3), lat_v.reshape(n, 3)], axis=-1
+                [
+                    geographic[:, 0].reshape(n, 3),
+                    geographic[:, 1].reshape(n, 3),
+                ],
+                axis=-1,
             )  # (N, 3, 2)
 
         with open(fname, "w") as fh:
diff --git a/src/geodef/geomap.py b/src/geodef/geomap.py
index 9b84d60..61e794d 100644
--- a/src/geodef/geomap.py
+++ b/src/geodef/geomap.py
@@ -15,8 +15,6 @@
 
 import numpy as np
 
-from geodef import transforms
-
 if TYPE_CHECKING:
     from geodef.data import GNSS
     from geodef.fault import Fault
@@ -109,15 +107,12 @@ def _patch_outlines_lonlat(fault: Fault) -> np.ndarray:
     verts = fault._vertices  # (N, 3, 3) as [e, n, u]
     assert verts is not None
     n_tri = verts.shape[0]
-    lon, lat, _ = transforms.enu2geod(
-        verts[:, :, 0].ravel(),
-        verts[:, :, 1].ravel(),
-        verts[:, :, 2].ravel(),
-        fault._ref_lat,
-        fault._ref_lon,
-        0.0,
+    geographic = fault.frame.to_geographic(
+        east=verts[:, :, 0].ravel(),
+        north=verts[:, :, 1].ravel(),
+        up=verts[:, :, 2].ravel(),
     )
-    lonlat = np.stack([lon, lat], axis=1).reshape(n_tri, 3, 2)
+    lonlat = geographic[:, :2].reshape(n_tri, 3, 2)
     return np.concatenate([lonlat, lonlat[:, :1, :]], axis=1)
 
 
diff --git a/src/geodef/geometry.py b/src/geodef/geometry.py
new file mode 100644
index 0000000..3e00fe9
--- /dev/null
+++ b/src/geodef/geometry.py
@@ -0,0 +1,349 @@
+"""Coordinate frames and geometry-array conversion functions."""
+
+from __future__ import annotations
+
+from collections.abc import Mapping
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.typing as npt
+
+from geodef import transforms
+from geodef.validation import check_finite_scalar, check_positive, check_range
+
+DEFAULT_PROJECTION = "wgs84-enu"
+"""Current geographic-to-local projection identifier."""
+
+PLANAR_PARAMETER_NAMES = (
+    "e0",
+    "n0",
+    "depth",
+    "strike",
+    "dip",
+    "length",
+    "width",
+)
+"""Order of the expert planar-geometry parameter vector."""
+
+__all__ = [
+    "DEFAULT_PROJECTION",
+    "PLANAR_PARAMETER_NAMES",
+    "LocalFrame",
+    "as_planar_vector",
+    "planar_parameter_dict",
+    "triangle_strike_dip",
+    "vertices_from_nodes",
+]
+
+
+@dataclass(frozen=True)
+class LocalFrame:
+    """A local East-North-Up frame tied to a geographic origin.
+
+    The only projection currently supported is ``"wgs84-enu"``: geographic
+    WGS84 coordinates are converted through Earth-centered Earth-fixed (ECEF)
+    coordinates and rotated into the tangent ENU frame at the origin.
+
+    Args:
+        origin_lat: Origin latitude in degrees.
+        origin_lon: Origin longitude in degrees.
+        origin_alt: Origin ellipsoidal altitude in meters.
+        projection: Geographic-to-local projection identifier. Currently only
+            ``"wgs84-enu"`` is supported.
+
+    Raises:
+        ValueError: If an origin value is non-finite or out of range, or the
+            projection is unsupported.
+    """
+
+    origin_lat: float
+    origin_lon: float
+    origin_alt: float = 0.0
+    projection: str = DEFAULT_PROJECTION
+
+    def __post_init__(self) -> None:
+        for name in ("origin_lat", "origin_lon", "origin_alt"):
+            value = float(getattr(self, name))
+            check_finite_scalar(
+                name, value, unit="degrees" if name != "origin_alt" else "meters"
+            )
+            object.__setattr__(self, name, value)
+        check_range("origin_lat", self.origin_lat, -90.0, 90.0, unit="degrees")
+        check_range("origin_lon", self.origin_lon, -360.0, 360.0, unit="degrees")
+        if self.projection != DEFAULT_PROJECTION:
+            raise ValueError(
+                f"projection must be {DEFAULT_PROJECTION!r}; got {self.projection!r}"
+            )
+
+    def to_enu(
+        self,
+        *,
+        lon: npt.ArrayLike,
+        lat: npt.ArrayLike,
+        alt: npt.ArrayLike,
+    ) -> np.ndarray:
+        """Convert geographic coordinates to ENU meters in this frame.
+
+        Args:
+            lon: Longitude in degrees.
+            lat: Latitude in degrees.
+            alt: Ellipsoidal altitude in meters.
+
+        Returns:
+            Array with final axis ``[east_m, north_m, up_m]``.
+        """
+        lon_array, lat_array, alt_array = np.broadcast_arrays(
+            np.asarray(lon, dtype=float),
+            np.asarray(lat, dtype=float),
+            np.asarray(alt, dtype=float),
+        )
+        if not np.all(
+            np.isfinite(np.stack([lon_array, lat_array, alt_array], axis=-1))
+        ):
+            raise ValueError("lon, lat, and alt must contain only finite values")
+        shape = lon_array.shape
+        east, north, up = transforms.geod2enu(
+            lat_array.ravel(),
+            lon_array.ravel(),
+            alt_array.ravel(),
+            self.origin_lat,
+            self.origin_lon,
+            self.origin_alt,
+        )
+        return np.stack([east, north, up], axis=-1).reshape((*shape, 3))
+
+    def to_geographic(
+        self,
+        *,
+        east: npt.ArrayLike,
+        north: npt.ArrayLike,
+        up: npt.ArrayLike,
+    ) -> np.ndarray:
+        """Convert ENU meters in this frame to geographic coordinates.
+
+        Args:
+            east: East coordinate in meters.
+            north: North coordinate in meters.
+            up: Up coordinate in meters.
+
+        Returns:
+            Array with final axis ``[longitude_degrees, latitude_degrees,
+            altitude_m]``.
+        """
+        east_array, north_array, up_array = np.broadcast_arrays(
+            np.asarray(east, dtype=float),
+            np.asarray(north, dtype=float),
+            np.asarray(up, dtype=float),
+        )
+        if not np.all(
+            np.isfinite(np.stack([east_array, north_array, up_array], axis=-1))
+        ):
+            raise ValueError("east, north, and up must contain only finite values")
+        shape = east_array.shape
+        lat, lon, alt = transforms.enu2geod(
+            east_array.ravel(),
+            north_array.ravel(),
+            up_array.ravel(),
+            self.origin_lat,
+            self.origin_lon,
+            self.origin_alt,
+        )
+        return np.stack([lon, lat, alt], axis=-1).reshape((*shape, 3))
+
+    def transform_enu(
+        self, coordinates: npt.ArrayLike, *, target: LocalFrame
+    ) -> np.ndarray:
+        """Explicitly re-express ENU coordinates in another local frame.
+
+        Args:
+            coordinates: ENU coordinates with shape ``(..., 3)``.
+            target: Destination frame.
+
+        Returns:
+            Coordinates in ``target`` with the same shape.
+
+        Raises:
+            ValueError: If ``coordinates`` does not have a final axis of 3.
+        """
+        array = np.asarray(coordinates, dtype=float)
+        if array.ndim == 0 or array.shape[-1] != 3:
+            raise ValueError(f"coordinates must have shape (..., 3), got {array.shape}")
+        geographic = self.to_geographic(
+            east=array[..., 0], north=array[..., 1], up=array[..., 2]
+        )
+        return target.to_enu(
+            lon=geographic[..., 0],
+            lat=geographic[..., 1],
+            alt=geographic[..., 2],
+        )
+
+    def is_compatible(self, other: LocalFrame) -> bool:
+        """Return whether ``other`` is exactly the same coordinate frame."""
+        return self == other
+
+    def require_compatible(self, other: LocalFrame) -> None:
+        """Raise if ``other`` is not exactly the same coordinate frame.
+
+        Args:
+            other: Frame that will be combined with this one.
+
+        Raises:
+            ValueError: If the origins or projection differ.
+        """
+        if not self.is_compatible(other):
+            raise ValueError(
+                "incompatible local frames: explicitly transform coordinates "
+                f"from {other!r} to {self!r} before combining them"
+            )
+
+
+def as_planar_vector(
+    parameters: npt.ArrayLike | Mapping[str, float],
+) -> np.ndarray:
+    """Validate and return the expert planar-geometry vector.
+
+    Args:
+        parameters: Either a mapping with the keys in
+            :data:`PLANAR_PARAMETER_NAMES` or an array already ordered as
+            ``[e0, n0, depth, strike, dip, length, width]``.
+
+    Returns:
+        Validated float array with shape ``(7,)``.
+
+    Raises:
+        ValueError: If keys, shape, finiteness, or physical ranges are invalid.
+    """
+    if isinstance(parameters, Mapping):
+        missing = [name for name in PLANAR_PARAMETER_NAMES if name not in parameters]
+        extra = [name for name in parameters if name not in PLANAR_PARAMETER_NAMES]
+        if missing or extra:
+            raise ValueError(
+                f"planar parameters have missing keys {missing} and extra keys {extra}"
+            )
+        values = np.array([parameters[name] for name in PLANAR_PARAMETER_NAMES])
+    else:
+        values = np.asarray(parameters, dtype=float)
+    if values.shape != (7,):
+        raise ValueError(f"planar parameters must have shape (7,), got {values.shape}")
+    if not np.all(np.isfinite(values)):
+        bad_index = int(np.flatnonzero(~np.isfinite(values))[0])
+        name = PLANAR_PARAMETER_NAMES[bad_index]
+        raise ValueError(f"{name} must be finite")
+
+    e0, n0, depth, strike, dip, length, width = map(float, values)
+    check_finite_scalar("e0", e0, unit="meters")
+    check_finite_scalar("n0", n0, unit="meters")
+    if depth < 0.0:
+        raise ValueError(f"depth must be non-negative (meters); got {depth:g}")
+    if not 0.0 <= strike < 360.0:
+        raise ValueError(f"strike must lie in [0, 360) degrees; got {strike:g}")
+    check_range("dip", dip, 0.0, 90.0, unit="degrees")
+    check_positive("length", length, unit="meters")
+    check_positive("width", width, unit="meters")
+    return np.array(values, dtype=float, copy=True)
+
+
+def planar_parameter_dict(
+    parameters: npt.ArrayLike | Mapping[str, float],
+) -> dict[str, float]:
+    """Return named planar parameters from a mapping or expert vector.
+
+    Args:
+        parameters: Planar parameter mapping or seven-element expert vector.
+
+    Returns:
+        Dictionary keyed by :data:`PLANAR_PARAMETER_NAMES`.
+    """
+    values = as_planar_vector(parameters)
+    return dict(zip(PLANAR_PARAMETER_NAMES, map(float, values), strict=True))
+
+
+def vertices_from_nodes(
+    nodes_enu: npt.ArrayLike, triangles: npt.ArrayLike
+) -> np.ndarray:
+    """Expand shared ENU nodes into per-triangle vertices.
+
+    Args:
+        nodes_enu: Shared nodes with shape ``(M, 3)``.
+        triangles: Integer node indices with shape ``(N, 3)``.
+
+    Returns:
+        Per-triangle ENU vertices with shape ``(N, 3, 3)``.
+
+    Raises:
+        ValueError: If array shapes, finiteness, or indices are invalid.
+    """
+    nodes = np.asarray(nodes_enu, dtype=float)
+    connectivity = np.asarray(triangles, dtype=int)
+    if nodes.ndim != 2 or nodes.shape[1] != 3:
+        raise ValueError(
+            f"nodes_enu node array must have shape (M, 3), got {nodes.shape}"
+        )
+    if not np.all(np.isfinite(nodes)):
+        raise ValueError("nodes_enu must contain only finite values")
+    if connectivity.ndim != 2 or connectivity.shape[1] != 3:
+        raise ValueError(f"triangles must have shape (N, 3), got {connectivity.shape}")
+    if connectivity.size and (
+        connectivity.min() < 0 or connectivity.max() >= nodes.shape[0]
+    ):
+        raise ValueError("triangles index out of range for nodes_enu")
+    return np.array(nodes[connectivity], copy=True)
+
+
+def triangle_strike_dip(
+    vertices_enu: npt.ArrayLike,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Derive right-hand-rule strike and dip from triangle vertices.
+
+    Args:
+        vertices_enu: Per-triangle vertices with shape ``(N, 3, 3)``.
+
+    Returns:
+        ``(strike_degrees, dip_degrees)`` arrays, each shape ``(N,)``.
+
+    Raises:
+        ValueError: If vertices have the wrong shape or contain non-finite
+            values.
+    """
+    vertices = np.asarray(vertices_enu, dtype=float)
+    if vertices.ndim != 3 or vertices.shape[1:] != (3, 3):
+        raise ValueError(
+            f"vertices_enu must have shape (N, 3, 3), got {vertices.shape}"
+        )
+    if not np.all(np.isfinite(vertices)):
+        raise ValueError("vertices_enu must contain only finite values")
+    edge1 = vertices[:, 1, :] - vertices[:, 0, :]
+    edge2 = vertices[:, 2, :] - vertices[:, 0, :]
+    normal = np.cross(edge1, edge2)
+    normal[normal[:, 2] < 0.0] *= -1.0
+    normal /= np.maximum(np.linalg.norm(normal, axis=1, keepdims=True), 1e-30)
+
+    dip = np.degrees(np.arccos(np.clip(np.abs(normal[:, 2]), 0.0, 1.0)))
+    strike = np.zeros(vertices.shape[0], dtype=float)
+    dipping = dip > 0.1
+    updip_azimuth = (
+        np.degrees(np.arctan2(normal[dipping, 0], normal[dipping, 1])) % 360.0
+    )
+    strike[dipping] = (updip_azimuth + 90.0) % 360.0
+    return strike, dip
+
+
+def _resolve_frame(
+    frame: LocalFrame | None,
+    ref_lat: float | None,
+    ref_lon: float | None,
+) -> LocalFrame:
+    """Resolve an explicit or legacy local frame."""
+    if (ref_lat is None) != (ref_lon is None):
+        raise ValueError("ref_lat and ref_lon must be provided together")
+    legacy = (
+        LocalFrame(ref_lat, ref_lon)
+        if ref_lat is not None and ref_lon is not None
+        else None
+    )
+    if frame is not None and legacy is not None:
+        frame.require_compatible(legacy)
+    selected = frame if frame is not None else legacy
+    if selected is None:
+        raise ValueError("provide frame or both ref_lat and ref_lon")
+    return selected
diff --git a/src/geodef/greens.py b/src/geodef/greens.py
index 2ae9944..fd4d681 100644
--- a/src/geodef/greens.py
+++ b/src/geodef/greens.py
@@ -3,7 +3,7 @@
 Combines displacement/strain Green's matrix construction (from okada_greens)
 with patch grid generation, component Green's functions, and Laplacian
 regularization operators (from okada_utils). Also provides the polymorphic
-``greens()`` function for assembling projected Green's matrices from
+``matrix()`` function for assembling projected Green's matrices from
 ``Fault`` and ``DataSet`` objects.
 """
 
@@ -19,6 +19,7 @@
 
 from geodef import backend, okada85, transforms, tri
 from geodef import cache as _cache
+from geodef.geometry import LocalFrame
 
 if TYPE_CHECKING:
     from geodef.data import DataSet
@@ -491,6 +492,8 @@ def tri_displacement_greens(
     depth: np.ndarray,
     vertices: np.ndarray,
     nu: float = 0.25,
+    *,
+    frame: LocalFrame | None = None,
 ) -> np.ndarray:
     """Build displacement Green's matrix for triangular fault patches.
 
@@ -506,6 +509,8 @@ def tri_displacement_greens(
         depth: Patch centroid depths (npatch,), positive down.
         vertices: Triangle vertices in local ENU, shape (npatch, 3, 3).
         nu: Poisson's ratio.
+        frame: Local frame defining ``vertices``. Defaults to the legacy
+            mean-centroid frame inferred from ``lat0`` and ``lon0``.
 
     Returns:
         G matrix of shape (3*nobs, 2*npatch). Columns ``[:npatch]`` are
@@ -516,12 +521,13 @@ def tri_displacement_greens(
     alt = np.zeros_like(lon)
     nobs = lat.shape[0]
     npatch = vertices.shape[0]
-    ref_lat = float(np.mean(lat0))
-    ref_lon = float(np.mean(lon0))
-
-    # Convert observation points to local ENU relative to fault centroid
-    obs_e, obs_n, _ = transforms.geod2enu(lat, lon, alt, ref_lat, ref_lon, 0.0)
-    obs = np.column_stack([obs_e, obs_n, np.zeros(nobs)])
+    selected_frame = (
+        LocalFrame(float(np.mean(lat0)), float(np.mean(lon0)))
+        if frame is None
+        else frame
+    )
+    obs_enu = selected_frame.to_enu(lon=lon, lat=lat, alt=alt)
+    obs = np.column_stack([obs_enu[:, 0], obs_enu[:, 1], np.zeros(nobs)])
 
     G = np.zeros((3 * nobs, 2 * npatch))
 
@@ -556,6 +562,8 @@ def tri_strain_greens(
     vertices: np.ndarray,
     nu: float = 0.25,
     obs_depth: np.ndarray | None = None,
+    *,
+    frame: LocalFrame | None = None,
 ) -> np.ndarray:
     """Build strain Green's matrix for triangular fault patches.
 
@@ -575,6 +583,8 @@ def tri_strain_greens(
             observations are at the surface. If provided, the z-coordinate
             of observation points is set to ``-obs_depth`` (negative = below
             surface in the ENU frame used by TDstrainHS).
+        frame: Local frame defining ``vertices``. Defaults to the legacy
+            mean-centroid frame inferred from ``lat0`` and ``lon0``.
 
     Returns:
         G matrix of shape (6*nobs, 2*npatch). Columns ``[:npatch]`` are
@@ -585,15 +595,17 @@ def tri_strain_greens(
     alt = np.zeros_like(lon)
     nobs = lat.shape[0]
     npatch = vertices.shape[0]
-    ref_lat = float(np.mean(lat0))
-    ref_lon = float(np.mean(lon0))
-
-    obs_e, obs_n, _ = transforms.geod2enu(lat, lon, alt, ref_lat, ref_lon, 0.0)
+    selected_frame = (
+        LocalFrame(float(np.mean(lat0)), float(np.mean(lon0)))
+        if frame is None
+        else frame
+    )
+    obs_enu = selected_frame.to_enu(lon=lon, lat=lat, alt=alt)
     if obs_depth is not None:
         obs_z = -np.asarray(obs_depth, dtype=float)
     else:
         obs_z = np.zeros(nobs)
-    obs = np.column_stack([obs_e, obs_n, obs_z])
+    obs = np.column_stack([obs_enu[:, 0], obs_enu[:, 1], obs_z])
 
     G = np.zeros((6 * nobs, 2 * npatch))
 
@@ -641,6 +653,10 @@ def _build_greens_key(fault: Fault, data: DataSet) -> dict:
         key["fault_width"] = fault._width
     if fault._vertices is not None:
         key["fault_vertices"] = fault._vertices
+        key["frame_origin_lat"] = fault.frame.origin_lat
+        key["frame_origin_lon"] = fault.frame.origin_lon
+        key["frame_origin_alt"] = fault.frame.origin_alt
+        key["frame_projection"] = fault.frame.projection
     if isinstance(data, InSAR):
         key["look_e"] = data._look_e
         key["look_n"] = data._look_n
@@ -657,6 +673,7 @@ def select_slip_columns(
     rake: float | None = None,
     fault_strike: np.ndarray | None = None,
     slip_azimuth: float | None = None,
+    plate_rake: float | np.ndarray | None = None,
 ) -> np.ndarray:
     """Project a two-component Green's matrix onto the requested slip basis.
 
@@ -669,15 +686,19 @@ def select_slip_columns(
         n_patches: Number of fault patches N.
         components: ``'both'`` (no reduction), ``'strike'``, ``'dip'``,
             ``'rake'`` (fixed rake, all patches), or ``'azimuth'`` (fixed
-            geographic slip azimuth, per-patch local rake).
+            geographic slip azimuth, per-patch local rake), or ``'plate'``
+            (two components parallel/perpendicular to a per-patch plate rake).
         rake: Fixed rake angle in degrees, required for ``'rake'``.
         fault_strike: Per-patch strike angles in degrees, shape (N,),
             required for ``'azimuth'``.
         slip_azimuth: Geographic slip azimuth in degrees CW from North,
             required for ``'azimuth'``.
+        plate_rake: Large-scale plate direction expressed as local rake in
+            degrees, scalar or shape (N,), required for ``'plate'``.
 
     Returns:
-        Reduced matrix: shape (M, 2*N) for ``'both'``, else (M, N).
+        Reduced matrix: shape (M, 2*N) for ``'both'`` and ``'plate'``, else
+        (M, N).
 
     Raises:
         ValueError: If required angles for the chosen basis are missing.
@@ -688,26 +709,40 @@ def select_slip_columns(
         return G_full[:, :n_patches]
     if components == "dip":
         return G_full[:, n_patches:]
+    if components == "plate":
+        if plate_rake is None:
+            raise ValueError("components='plate' requires plate_rake")
+        theta = np.deg2rad(np.broadcast_to(plate_rake, (n_patches,)))
+        cosine = np.cos(theta)
+        sine = np.sin(theta)
+        strike = G_full[:, :n_patches]
+        dip = G_full[:, n_patches:]
+        parallel = strike * cosine + dip * sine
+        perpendicular = -strike * sine + dip * cosine
+        return np.hstack([parallel, perpendicular])
     if components == "rake":
         if rake is None:
             raise ValueError("components='rake' requires a rake angle in degrees")
         theta = np.deg2rad(rake)  # scalar
-    else:  # azimuth: per-patch local rake = slip_azimuth - strike_i
+    elif components == "azimuth":
         if fault_strike is None or slip_azimuth is None:
             raise ValueError(
                 "components='azimuth' requires fault_strike and slip_azimuth"
             )
         theta = np.deg2rad(slip_azimuth - fault_strike)  # shape (N,)
+    else:
+        raise ValueError(f"Unknown slip components {components!r}")
     return G_full[:, :n_patches] * np.cos(theta) + G_full[:, n_patches:] * np.sin(theta)
 
 
-def greens(
+def matrix(
     fault: Fault,
     datasets: DataSet | list[DataSet],
     *,
     components: str = "both",
     rake: float | None = None,
     slip_azimuth: float | None = None,
+    plate_rake: float | np.ndarray | None = None,
 ) -> np.ndarray:
     """Build a projected Green's matrix for one or more datasets.
 
@@ -720,11 +755,15 @@ def greens(
         datasets: A single ``DataSet`` or a list of them.
         components: Slip basis for the returned columns: ``'both'``
             (default, ``2*N`` columns), ``'strike'``, ``'dip'``, ``'rake'``
-            (fixed rake), or ``'azimuth'`` (fixed geographic slip azimuth).
-            The reduction uses the same semantics as :func:`geodef.invert`.
+            (fixed rake), ``'azimuth'`` (fixed geographic slip azimuth), or
+            ``'plate'`` (plate-rake-parallel/perpendicular coordinates).
+            The reduction uses the same semantics as
+            :func:`geodef.invert.solve`.
         rake: Fixed rake angle in degrees, required for ``components='rake'``.
         slip_azimuth: Geographic slip azimuth in degrees CW from North,
             required for ``components='azimuth'``.
+        plate_rake: Large-scale direction as local rake, scalar or shape (N,),
+            required for ``components='plate'``.
 
     Returns:
         Projected Green's matrix. For a single dataset with M observations
@@ -742,7 +781,7 @@ def greens(
         key = _build_greens_key(fault, data)
         G_proj = _cache.cached_compute(
             key,
-            lambda: _project_greens(
+            lambda: project(
                 data, fault.greens_matrix(data.lat, data.lon, kind=data.greens_type)
             ),
         )
@@ -756,6 +795,7 @@ def greens(
         rake,
         fault_strike=fault.strike,
         slip_azimuth=slip_azimuth,
+        plate_rake=plate_rake,
     )
 
 
@@ -800,7 +840,7 @@ def stack_weights(datasets: DataSet | list[DataSet]) -> np.ndarray:
     return scipy.linalg.block_diag(*blocks)
 
 
-def _project_greens(data: DataSet, G_raw: np.ndarray) -> np.ndarray:
+def project(data: DataSet, G_raw: np.ndarray) -> np.ndarray:
     """Project a raw Green's matrix through a dataset's projection.
 
     For displacement Green's matrices (3 components per station), reshapes
@@ -837,6 +877,18 @@ def _project_greens(data: DataSet, G_raw: np.ndarray) -> np.ndarray:
     return G_proj
 
 
+def laplacian(fault: Fault) -> np.ndarray:
+    """Return the fault's patch Laplacian regularization matrix.
+
+    Args:
+        fault: Rectangular or triangular fault.
+
+    Returns:
+        Dense Laplacian matrix with shape ``(N, N)``.
+    """
+    return fault.laplacian
+
+
 def resolution(G: np.ndarray) -> np.ndarray:
     """Compute resolution matrix R = pinv(G) @ G.
 
diff --git a/src/geodef/invert.py b/src/geodef/invert.py
index a93ffcb..a6d1fbf 100644
--- a/src/geodef/invert.py
+++ b/src/geodef/invert.py
@@ -8,6 +8,7 @@
 
 import dataclasses
 import functools
+from collections.abc import Mapping
 from pathlib import Path
 from typing import TYPE_CHECKING
 
@@ -21,12 +22,19 @@
 from geodef import backend
 from geodef.data import DataSet
 from geodef.fault import Fault, moment_to_magnitude
-from geodef.greens import greens, select_slip_columns, stack_obs, stack_weights
+from geodef.geometry import (
+    LocalFrame,
+    _resolve_frame,
+    as_planar_vector,
+)
+from geodef.greens import matrix, select_slip_columns, stack_obs, stack_weights
+from geodef.slip import from_plate, from_rake, magnitude, unpack
+from geodef.slip import rake as slip_rake
 
 _VALID_METHODS = {"wls", "nnls", "bounded_ls", "constrained"}
 _VALID_SMOOTHING_STRINGS = {"laplacian", "damping", "stresskernel"}
 _VALID_STRENGTH_STRINGS = {"abic", "cv"}
-_VALID_COMPONENTS = {"both", "strike", "dip", "rake", "azimuth"}
+_VALID_COMPONENTS = {"both", "strike", "dip", "rake", "azimuth", "plate"}
 
 # A bound may be a scalar (all parameters), an array of length n_components
 # (one value per slip component, broadcast over patches), or an array of
@@ -66,6 +74,9 @@ class InversionResult:
             North when ``components='azimuth'``, else ``None``. Each
             patch's effective local rake is ``slip_azimuth - strike_i``,
             so this correctly handles faults with varying strike.
+        plate_rake: Large-scale plate direction expressed as local rake per
+            patch when ``components='plate'``. The two solution blocks are
+            rake-parallel and rake-perpendicular.
     """
 
     slip: np.ndarray
@@ -81,6 +92,61 @@ class InversionResult:
     components: str
     rake: float | None = None
     slip_azimuth: float | None = None
+    plate_rake: np.ndarray | None = None
+    local_rake: np.ndarray | None = None
+
+    @property
+    def n_patches(self) -> int:
+        """Number of fault patches represented by the result."""
+        divisor = 2 if self.components in {"both", "plate"} else 1
+        return self.slip_vector.size // divisor
+
+    @property
+    def strike_slip(self) -> np.ndarray:
+        """Physical strike-slip component per patch."""
+        return self._physical_components()[0]
+
+    @property
+    def dip_slip(self) -> np.ndarray:
+        """Physical dip-slip component per patch."""
+        return self._physical_components()[1]
+
+    @property
+    def slip_magnitude(self) -> np.ndarray:
+        """Unsigned physical slip magnitude per patch."""
+        return magnitude(self.strike_slip, self.dip_slip)
+
+    @property
+    def slip_rake(self) -> np.ndarray:
+        """Physical local rake in degrees per patch."""
+        return slip_rake(self.strike_slip, self.dip_slip)
+
+    @property
+    def rake_parallel(self) -> np.ndarray:
+        """Plate-rake-parallel solution component per patch."""
+        if self.components != "plate":
+            raise AttributeError("rake_parallel requires components='plate'")
+        return self.slip_vector[: self.n_patches]
+
+    @property
+    def rake_perpendicular(self) -> np.ndarray:
+        """Plate-rake-perpendicular solution component per patch."""
+        if self.components != "plate":
+            raise AttributeError("rake_perpendicular requires components='plate'")
+        return self.slip_vector[self.n_patches :]
+
+    def _physical_components(self) -> tuple[np.ndarray, np.ndarray]:
+        """Convert the solved basis to physical strike/dip components."""
+        angle: float | np.ndarray | None
+        if self.components == "rake":
+            angle = self.rake
+        elif self.components == "azimuth":
+            angle = self.local_rake
+        elif self.components == "plate":
+            angle = self.plate_rake
+        else:
+            angle = None
+        return _physical_components(self.slip_vector, self.components, angle)
 
     # ------------------------------------------------------------------
     # I/O
@@ -122,6 +188,16 @@ def save(self, fname: str | Path) -> None:
             if self.slip_azimuth is None
             else np.array([self.slip_azimuth])
         )
+        plate_rake_arr = (
+            np.array([], dtype=float)
+            if self.plate_rake is None
+            else np.asarray(self.plate_rake, dtype=float)
+        )
+        local_rake_arr = (
+            np.array([], dtype=float)
+            if self.local_rake is None
+            else np.asarray(self.local_rake, dtype=float)
+        )
 
         arrays: dict = {
             "slip": self.slip,
@@ -137,6 +213,8 @@ def save(self, fname: str | Path) -> None:
             "components": np.array([self.components]),
             "rake": rake_arr,
             "slip_azimuth": slip_azimuth_arr,
+            "plate_rake": plate_rake_arr,
+            "local_rake": local_rake_arr,
         }
         if smoothing_arr is not None:
             arrays["smoothing_arr"] = smoothing_arr
@@ -173,6 +251,10 @@ def load(cls, fname: str | Path) -> "InversionResult":
             float(data["slip_azimuth"][0]) if "slip_azimuth" in data else float("nan")
         )
         slip_azimuth: float | None = None if np.isnan(raw_az) else raw_az
+        plate_rake = data["plate_rake"] if "plate_rake" in data else np.array([])
+        plate_rake_value = None if plate_rake.size == 0 else plate_rake
+        local_rake = data["local_rake"] if "local_rake" in data else np.array([])
+        local_rake_value = None if local_rake.size == 0 else local_rake
 
         return cls(
             slip=data["slip"],
@@ -188,6 +270,8 @@ def load(cls, fname: str | Path) -> "InversionResult":
             components=str(data["components"][0]),
             rake=rake,
             slip_azimuth=slip_azimuth,
+            plate_rake=plate_rake_value,
+            local_rake=local_rake_value,
         )
 
     def save_table(self, fname: str | Path, fault: "Fault") -> None:
@@ -263,8 +347,10 @@ def save_table(self, fname: str | Path, fault: "Fault") -> None:
             slip_col_names = ["slip_strike_m"]
         elif self.components == "dip":
             slip_col_names = ["slip_dip_m"]
-        else:  # rake or azimuth
+        elif self.components in {"rake", "azimuth"}:
             slip_col_names = ["slip_amplitude_m"]
+        else:  # plate
+            slip_col_names = ["slip_rake_parallel_m", "slip_rake_perpendicular_m"]
         slip_cols = "  ".join(slip_col_names)
         header_lines.append(f"{col_names}  {slip_cols}")
 
@@ -374,9 +460,11 @@ class GeometrySearchResult:
     """Result of a gradient-based nonlinear geometry search.
 
     Attributes:
+        fault: Optimal fault geometry.
+        frame: Local frame defining ``theta``.
         theta: Optimal geometry, full 7-vector
             ``[e0, n0, depth, strike, dip, length, width]`` in the local
-            Cartesian frame anchored at (ref_lat, ref_lon).
+            Cartesian :attr:`geometry.frame`.
         free: Names of the parameters that were optimized.
         slip: Slip solved linearly at the optimal geometry (inner solve).
         chi2: Weighted misfit ``r^T W r`` at the optimum.
@@ -388,6 +476,8 @@ class GeometrySearchResult:
         n_iterations: Number of optimizer iterations.
     """
 
+    fault: Fault
+    frame: LocalFrame
     theta: np.ndarray
     free: list[str]
     slip: np.ndarray
@@ -496,7 +586,12 @@ class LinearSystem:
         smoothing: Regularization type — ``'laplacian'``, ``'damping'``,
             ``'stresskernel'``, a custom matrix, or ``None``.
         components: Slip components to solve for: ``'both'`` (default),
-            ``'strike'``, or ``'dip'``.
+            ``'strike'``, ``'dip'``, ``'rake'``, ``'azimuth'``, or ``'plate'``.
+        rake: Constant local rake for ``components='rake'``.
+        slip_azimuth: Constant geographic direction for
+            ``components='azimuth'``.
+        plate_rake: Scalar or per-patch large-scale direction in local rake
+            coordinates for ``components='plate'``.
 
     Examples:
         >>> sys = LinearSystem(fault, [gnss, insar], smoothing='laplacian')
@@ -513,6 +608,7 @@ def __init__(
         components: str = "both",
         rake: float | None = None,
         slip_azimuth: float | None = None,
+        plate_rake: float | np.ndarray | None = None,
     ) -> None:
         if isinstance(datasets, DataSet):
             datasets = [datasets]
@@ -539,6 +635,13 @@ def __init__(
                 f"slip_azimuth is only used with components='azimuth', "
                 f"got components={components!r}"
             )
+        if components == "plate" and plate_rake is None:
+            raise ValueError("components='plate' requires plate_rake")
+        if plate_rake is not None and components != "plate":
+            raise ValueError(
+                "plate_rake is only used with components='plate', "
+                f"got components={components!r}"
+            )
 
         self.fault = fault
         self.datasets = datasets
@@ -546,13 +649,20 @@ def __init__(
         self.components = components
         self.rake = rake
         self.slip_azimuth = slip_azimuth
+        self.plate_rake = (
+            None
+            if plate_rake is None
+            else np.broadcast_to(
+                np.asarray(plate_rake, dtype=float), (fault.n_patches,)
+            ).copy()
+        )
 
         n_patches = fault.n_patches
-        n_components = 2 if components == "both" else 1
+        n_components = 2 if components in {"both", "plate"} else 1
         self._n_patches = n_patches
         self._n_params = n_components * n_patches
 
-        G_full = greens(fault, datasets)
+        G_full = matrix(fault, datasets)
         self.d = stack_obs(datasets)
         self.W = stack_weights(datasets)
         self.G = select_slip_columns(
@@ -562,6 +672,7 @@ def __init__(
             rake,
             fault_strike=fault.strike,
             slip_azimuth=slip_azimuth,
+            plate_rake=self.plate_rake,
         )
         self.G_w, self.d_w = _apply_weights(self.G, self.d, self.W)
         self.L: np.ndarray | None = (
@@ -573,6 +684,7 @@ def __init__(
                 components,
                 rake,
                 slip_azimuth,
+                self.plate_rake,
             )
             if smoothing is not None
             else None
@@ -818,7 +930,7 @@ def invert(
             bounds: Per-component slip bounds ``(lower, upper)``.
             method: Solver — ``'wls'``, ``'nnls'``, ``'bounded_ls'``,
                 or ``'constrained'``. Auto-selected from bounds if None.
-            smoothing_target: Reference model, shape (n_params,).
+            smoothing_target: Reference vector, shape ``(n_params,)``.
                 Regularizes toward this target instead of zero.
             constraints: Inequality constraints ``(C, d_ineq)`` such
                 that ``C @ m <= d_ineq``.
@@ -841,6 +953,7 @@ def invert(
             self._n_params,
             self.rake,
             self.slip_azimuth,
+            self.plate_rake,
         )
 
         exp_bounds = _expand_bounds(
@@ -880,13 +993,22 @@ def invert(
         reduced_chi2 = _compute_reduced_chi2(residuals, self.W, self._n_params)
         rms = float(np.sqrt(np.mean(residuals**2)))
 
-        if self.components == "both":
+        if self.components in {"both", "plate"}:
             slip = np.column_stack([m[: self._n_patches], m[self._n_patches :]])
-            slip_mag = np.sqrt(slip[:, 0] ** 2 + slip[:, 1] ** 2)
         else:
             slip = m.reshape(-1, 1)
-            slip_mag = np.abs(m)
-        moment = self.fault.moment(slip_mag)
+        basis_angle: float | np.ndarray | None
+        if self.components == "rake":
+            basis_angle = self.rake
+        elif self.components == "azimuth":
+            assert self.slip_azimuth is not None
+            basis_angle = self.slip_azimuth - self.fault.strike
+        elif self.components == "plate":
+            basis_angle = self.plate_rake
+        else:
+            basis_angle = None
+        strike_slip, dip_slip = _physical_components(m, self.components, basis_angle)
+        moment = self.fault.moment(magnitude(strike_slip, dip_slip))
         mw = moment_to_magnitude(moment)
 
         return InversionResult(
@@ -903,6 +1025,12 @@ def invert(
             components=self.components,
             rake=self.rake,
             slip_azimuth=self.slip_azimuth,
+            plate_rake=self.plate_rake,
+            local_rake=(
+                self.slip_azimuth - self.fault.strike
+                if self.slip_azimuth is not None
+                else None
+            ),
         )
 
     def lcurve(
@@ -1137,7 +1265,7 @@ def model_uncertainty(
 # ======================================================================
 
 
-def invert(
+def solve(
     fault: Fault,
     datasets: DataSet | list[DataSet],
     smoothing: str | np.ndarray | None = None,
@@ -1150,6 +1278,7 @@ def invert(
     slip_azimuth: float | None = None,
     constraints: tuple[np.ndarray, np.ndarray] | None = None,
     cv_folds: int = 5,
+    plate_rake: float | np.ndarray | None = None,
 ) -> InversionResult:
     """Invert geodetic data for fault slip.
 
@@ -1179,6 +1308,9 @@ def invert(
             required when ``components='azimuth'``. Each patch's
             effective local rake is ``slip_azimuth - strike_i``,
             so this correctly handles faults with varying strike.
+        plate_rake: Large-scale direction as a local rake angle, scalar or
+            shape (N,), required when ``components='plate'``. The solved
+            blocks are rake-parallel and rake-perpendicular.
         constraints: Inequality constraints ``(C, d_ineq)`` such that
             ``C @ m <= d_ineq``. Only used with ``method='constrained'``.
         cv_folds: Number of folds for cross-validation (default 5).
@@ -1189,7 +1321,15 @@ def invert(
     Raises:
         ValueError: For invalid arguments.
     """
-    sys = LinearSystem(fault, datasets, smoothing, components, rake, slip_azimuth)
+    sys = LinearSystem(
+        fault,
+        datasets,
+        smoothing,
+        components,
+        rake,
+        slip_azimuth,
+        plate_rake,
+    )
     return sys.invert(
         smoothing_strength, bounds, method, smoothing_target, constraints, cv_folds
     )
@@ -1251,6 +1391,7 @@ def lcurve(
     components: str = "both",
     rake: float | None = None,
     slip_azimuth: float | None = None,
+    plate_rake: float | np.ndarray | None = None,
 ) -> LCurveResult:
     """Sweep smoothing strength and compute the L-curve.
 
@@ -1267,11 +1408,15 @@ def lcurve(
             ``components='rake'``.
         slip_azimuth: Geographic slip azimuth in degrees, required when
             ``components='azimuth'``.
+        plate_rake: Local plate-rake direction, required when
+            ``components='plate'``.
 
     Returns:
         LCurveResult with sweep arrays and optimal lambda.
     """
-    sys = LinearSystem(fault, datasets, smoothing, components, rake, slip_azimuth)
+    sys = LinearSystem(
+        fault, datasets, smoothing, components, rake, slip_azimuth, plate_rake
+    )
     return sys.lcurve(smoothing_range, n, bounds, method)
 
 
@@ -1284,6 +1429,7 @@ def abic_curve(
     components: str = "both",
     rake: float | None = None,
     slip_azimuth: float | None = None,
+    plate_rake: float | np.ndarray | None = None,
 ) -> ABICCurveResult:
     """Sweep smoothing strength and compute the ABIC at each value.
 
@@ -1301,11 +1447,15 @@ def abic_curve(
             ``components='rake'``.
         slip_azimuth: Geographic slip azimuth in degrees, required when
             ``components='azimuth'``.
+        plate_rake: Local plate-rake direction, required when
+            ``components='plate'``.
 
     Returns:
         ABICCurveResult with sweep arrays and optimal lambda.
     """
-    sys = LinearSystem(fault, datasets, smoothing, components, rake, slip_azimuth)
+    sys = LinearSystem(
+        fault, datasets, smoothing, components, rake, slip_azimuth, plate_rake
+    )
     return sys.abic_curve(smoothing_range, n)
 
 
@@ -1410,12 +1560,35 @@ def _vp_kernel():
     return _vp_jitted["kernel"]
 
 
+def _fault_from_planar_vector(
+    theta: np.ndarray,
+    frame: LocalFrame,
+    n_length: int,
+    n_width: int,
+) -> Fault:
+    """Construct a planar fault from the local expert parameter vector."""
+    geographic = frame.to_geographic(east=theta[0], north=theta[1], up=0.0)
+    return Fault.planar(
+        lat=float(geographic[1]),
+        lon=float(geographic[0]),
+        depth=float(theta[2]),
+        strike=float(theta[3]),
+        dip=float(theta[4]),
+        length=float(theta[5]),
+        width=float(theta[6]),
+        n_length=n_length,
+        n_width=n_width,
+        frame=frame,
+    )
+
+
 def geometry_search(
-    theta0: np.ndarray,
+    theta0: np.ndarray | Mapping[str, float],
     datasets: DataSet | list[DataSet],
     *,
-    ref_lat: float,
-    ref_lon: float,
+    ref_lat: float | None = None,
+    ref_lon: float | None = None,
+    frame: LocalFrame | None = None,
     free: list[str] | None = None,
     bounds: dict[str, tuple[float, float]] | None = None,
     n_length: int = 1,
@@ -1439,13 +1612,15 @@ def geometry_search(
     Requires the JAX backend (``geodef.backend.set_backend('jax')``).
 
     Args:
-        theta0: Starting geometry ``[e0, n0, depth, strike, dip, length,
-            width]``; ``e0``/``n0`` are centroid offsets in meters from
-            (ref_lat, ref_lon).
+        theta0: Starting parameter mapping, or expert array
+            ``[east, north, depth, strike, dip, length, width]``. Requires
+            ``frame`` or ``ref_lat``/``ref_lon``.
         datasets: One or more displacement datasets (GNSS, InSAR,
             Vertical).
         ref_lat: Latitude anchoring the local Cartesian frame.
         ref_lon: Longitude anchoring the local Cartesian frame.
+        frame: Explicit local frame for array ``theta0``. Mutually exclusive
+            with an incompatible legacy ``ref_lat``/``ref_lon`` origin.
         free: Names of parameters to optimize (subset of ``e0, n0,
             depth, strike, dip, length, width``). Default: all seven.
         bounds: Optional per-parameter ``(lower, upper)`` bounds, keyed
@@ -1461,8 +1636,8 @@ def geometry_search(
         nu: Poisson's ratio.
 
     Returns:
-        GeometrySearchResult with the optimal geometry, inner slip,
-        misfit, and a Gauss-Newton covariance for the free parameters.
+        GeometrySearchResult with optimal ``fault``, expert ``theta``, frame,
+        inner slip, misfit, and a Gauss-Newton covariance.
 
     Raises:
         RuntimeError: If the JAX backend is not active.
@@ -1476,8 +1651,6 @@ def geometry_search(
         )
     import jax.numpy as jnp
 
-    from geodef import transforms
-
     if isinstance(datasets, DataSet):
         datasets = [datasets]
     if free is None:
@@ -1493,22 +1666,13 @@ def geometry_search(
             f"'dip', got {components!r}"
         )
 
-    theta0 = np.asarray(theta0, dtype=float)
+    frame = _resolve_frame(frame, ref_lat, ref_lon)
+    theta0 = as_planar_vector(theta0)
     free_idx = np.array([_THETA_NAMES.index(name) for name in free])
 
     # Template system provides the stacked data, weights, and (fixed)
     # regularization operator; its Green's matrix is not used.
-    template = Fault.planar(
-        lat=ref_lat,
-        lon=ref_lon,
-        depth=theta0[2],
-        strike=theta0[3],
-        dip=theta0[4],
-        length=theta0[5],
-        width=theta0[6],
-        n_length=n_length,
-        n_width=n_width,
-    )
+    template = _fault_from_planar_vector(theta0, frame, n_length, n_width)
     sys = LinearSystem(template, datasets, smoothing, components)
     n_patches = n_length * n_width
     col_start, col_stop = {
@@ -1519,11 +1683,13 @@ def geometry_search(
 
     e_parts, n_parts = [], []
     for ds in datasets:
-        e_ds, n_ds, _ = transforms.geod2enu(
-            ds.lat, ds.lon, np.zeros(ds.n_stations), ref_lat, ref_lon, 0.0
+        enu = frame.to_enu(
+            lon=ds.lon,
+            lat=ds.lat,
+            alt=np.full(ds.n_stations, frame.origin_alt),
         )
-        e_parts.append(e_ds)
-        n_parts.append(n_ds)
+        e_parts.append(enu[:, 0])
+        n_parts.append(enu[:, 1])
     e_obs = np.concatenate(e_parts)
     n_obs = np.concatenate(n_parts)
 
@@ -1583,8 +1749,11 @@ def scipy_objective(x: np.ndarray) -> tuple[float, np.ndarray]:
 
     theta_opt = theta0.copy()
     theta_opt[free_idx] = np.asarray(opt.x, dtype=float)
+    fault_opt = _fault_from_planar_vector(theta_opt, frame, n_length, n_width)
 
     return GeometrySearchResult(
+        fault=fault_opt,
+        frame=frame,
         theta=theta_opt,
         free=list(free),
         slip=backend.to_numpy(m),
@@ -1623,6 +1792,7 @@ def dataset_diagnostics(
         result.components,
         result.rake,
         result.slip_azimuth,
+        result.plate_rake,
     )
     return sys.dataset_diagnostics(result)
 
@@ -1661,6 +1831,7 @@ def model_covariance(
         result.components,
         result.rake,
         result.slip_azimuth,
+        result.plate_rake,
     )
     return sys.model_covariance(result, kind=kind)
 
@@ -1693,6 +1864,7 @@ def model_resolution(
         result.components,
         result.rake,
         result.slip_azimuth,
+        result.plate_rake,
     )
     return sys.model_resolution(result)
 
@@ -1724,6 +1896,7 @@ def model_uncertainty(
         result.components,
         result.rake,
         result.slip_azimuth,
+        result.plate_rake,
     )
     return sys.model_uncertainty(result, kind=kind)
 
@@ -1733,6 +1906,30 @@ def model_uncertainty(
 # ======================================================================
 
 
+def _physical_components(
+    vector: np.ndarray,
+    components: str,
+    basis_angle: float | np.ndarray | None,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Convert a solved basis vector to physical strike/dip components."""
+    if components == "both":
+        return unpack(vector)
+    if components == "strike":
+        return vector, np.zeros_like(vector)
+    if components == "dip":
+        return np.zeros_like(vector), vector
+    if components in {"rake", "azimuth"}:
+        if basis_angle is None:
+            raise ValueError(f"{components} result is missing angle metadata")
+        return from_rake(vector, basis_angle)
+    if components == "plate":
+        if basis_angle is None:
+            raise ValueError("plate result is missing plate_rake metadata")
+        parallel, perpendicular = unpack(vector)
+        return from_plate(parallel, perpendicular, basis_angle)
+    raise ValueError(f"Unknown slip components {components!r}")
+
+
 def _validate_args(
     datasets: list[DataSet],
     components: str,
@@ -1744,6 +1941,7 @@ def _validate_args(
     n_params: int,
     rake: float | None = None,
     slip_azimuth: float | None = None,
+    plate_rake: np.ndarray | None = None,
 ) -> None:
     """Validate invert() arguments."""
     for ds in datasets:
@@ -1771,6 +1969,13 @@ def _validate_args(
             f"slip_azimuth is only used with components='azimuth', "
             f"got components={components!r}"
         )
+    if components == "plate" and plate_rake is None:
+        raise ValueError("components='plate' requires plate_rake")
+    if plate_rake is not None and components != "plate":
+        raise ValueError(
+            "plate_rake is only used with components='plate', "
+            f"got components={components!r}"
+        )
 
     if method is not None and method not in _VALID_METHODS:
         raise ValueError(f"method must be one of {_VALID_METHODS}, got {method!r}")
@@ -1838,6 +2043,7 @@ def _build_smoothing_matrix(
     components: str,
     rake: float | None = None,
     slip_azimuth: float | None = None,
+    plate_rake: np.ndarray | None = None,
 ) -> np.ndarray:
     """Build the regularization matrix L.
 
@@ -1850,6 +2056,7 @@ def _build_smoothing_matrix(
         rake: Fixed rake angle, used when ``components='rake'``.
         slip_azimuth: Fixed geographic slip azimuth, used when
             ``components='azimuth'``.
+        plate_rake: Per-patch plate rake, used when ``components='plate'``.
 
     Returns:
         Regularization matrix with n_params columns.
@@ -1875,6 +2082,7 @@ def _build_smoothing_matrix(
             rake,
             fault_strike=fault.strike,
             slip_azimuth=slip_azimuth,
+            plate_rake=plate_rake,
         )
 
     raise ValueError(f"Unknown smoothing type: {smoothing!r}")
diff --git a/src/geodef/mesh.py b/src/geodef/mesh.py
index d758cf1..296e3bf 100644
--- a/src/geodef/mesh.py
+++ b/src/geodef/mesh.py
@@ -24,6 +24,7 @@
 from scipy import interpolate
 
 from geodef import transforms
+from geodef.geometry import LocalFrame
 
 logger = logging.getLogger(__name__)
 
@@ -42,18 +43,25 @@ class Mesh:
         lat: Node latitudes, shape (N,).
         depth: Node depths in meters (positive down), shape (N,).
         triangles: Triangle connectivity as indices into nodes, shape (M, 3).
+        frame: Local frame used by :meth:`vertices_enu`. Defaults to a
+            ``wgs84-enu`` frame at the mean node latitude and longitude.
     """
 
     lon: np.ndarray
     lat: np.ndarray
     depth: np.ndarray
     triangles: np.ndarray
+    frame: LocalFrame | None = None
 
     def __post_init__(self) -> None:
         object.__setattr__(self, "lon", np.asarray(self.lon, dtype=float))
         object.__setattr__(self, "lat", np.asarray(self.lat, dtype=float))
         object.__setattr__(self, "depth", np.asarray(self.depth, dtype=float))
         object.__setattr__(self, "triangles", np.asarray(self.triangles, dtype=int))
+        if self.frame is None:
+            origin_lat = float(np.mean(self.lat)) if self.lat.size else 0.0
+            origin_lon = float(np.mean(self.lon)) if self.lon.size else 0.0
+            object.__setattr__(self, "frame", LocalFrame(origin_lat, origin_lon))
 
     @property
     def n_nodes(self) -> int:
@@ -107,9 +115,7 @@ def validate(self) -> "ValidationReport":
                     f"{tiny.size} degenerate triangle(s) with area < 1 m^2 "
                     f"(first indices {tiny[:5].tolist()})",
                 )
-            ref_lat = float(np.mean(self.lat))
-            ref_lon = float(np.mean(self.lon))
-            verts = self.vertices_enu(ref_lat, ref_lon)
+            verts = self.vertices_enu()
             edges = np.stack(
                 [
                     np.linalg.norm(verts[:, 1] - verts[:, 0], axis=1),
@@ -142,42 +148,77 @@ def validate(self) -> "ValidationReport":
     @property
     def areas(self) -> np.ndarray:
         """Triangle areas in m^2, shape (M,)."""
-        ref_lat = float(np.mean(self.lat))
-        ref_lon = float(np.mean(self.lon))
-        verts = self.vertices_enu(ref_lat, ref_lon)
+        verts = self.vertices_enu()
         edge1 = verts[:, 1, :] - verts[:, 0, :]
         edge2 = verts[:, 2, :] - verts[:, 0, :]
         return 0.5 * np.linalg.norm(np.cross(edge1, edge2), axis=1)
 
-    def vertices_enu(self, ref_lat: float, ref_lon: float) -> np.ndarray:
+    def vertices_enu(
+        self,
+        ref_lat: float | None = None,
+        ref_lon: float | None = None,
+        *,
+        frame: LocalFrame | None = None,
+    ) -> np.ndarray:
         """Triangle vertices in local ENU meters, shape (M, 3, 3).
 
         Each triangle has 3 vertices, each with [east, north, up] coordinates
         relative to the reference point. Depth is converted to up (z = -depth).
 
         Args:
-            ref_lat: Reference latitude for ENU origin.
-            ref_lon: Reference longitude for ENU origin.
+            ref_lat: Legacy reference latitude for an alternate ENU origin.
+            ref_lon: Legacy reference longitude for an alternate ENU origin.
+            frame: Explicit alternate frame. Mutually exclusive with
+                ``ref_lat`` and ``ref_lon``. Defaults to :attr:`frame`.
 
         Returns:
             Array of shape (M, 3, 3) suitable for ``Fault.__init__(vertices=...)``.
+
+        Raises:
+            ValueError: If only one legacy origin coordinate is supplied, or
+                legacy origin coordinates and ``frame`` are mixed.
         """
-        e, n, u = transforms.geod2enu(
-            self.lat,
-            self.lon,
-            -self.depth,
-            ref_lat,
-            ref_lon,
-            0.0,
+        if frame is not None and (ref_lat is not None or ref_lon is not None):
+            raise ValueError("provide either frame or ref_lat/ref_lon, not both")
+        if (ref_lat is None) != (ref_lon is None):
+            raise ValueError("ref_lat and ref_lon must be provided together")
+        if frame is not None:
+            selected_frame = frame
+        elif ref_lat is not None and ref_lon is not None:
+            selected_frame = LocalFrame(ref_lat, ref_lon)
+        else:
+            assert self.frame is not None
+            selected_frame = self.frame
+        nodes = selected_frame.to_enu(
+            lon=self.lon,
+            lat=self.lat,
+            alt=-self.depth,
         )
         tri = self.triangles
         verts = np.empty((self.n_triangles, 3, 3), dtype=float)
         for k in range(3):
-            verts[:, k, 0] = e[tri[:, k]]
-            verts[:, k, 1] = n[tri[:, k]]
-            verts[:, k, 2] = u[tri[:, k]]
+            verts[:, k, :] = nodes[tri[:, k]]
         return verts
 
+    def to_frame(self, frame: LocalFrame) -> "Mesh":
+        """Return this geographic mesh with a different default local frame.
+
+        Args:
+            frame: Destination frame for :meth:`vertices_enu`.
+
+        Returns:
+            Mesh with unchanged geographic nodes and connectivity.
+        """
+        if self.frame == frame:
+            return self
+        return Mesh(
+            lon=self.lon,
+            lat=self.lat,
+            depth=self.depth,
+            triangles=self.triangles,
+            frame=frame,
+        )
+
     # ------------------------------------------------------------------
     # I/O
     # ------------------------------------------------------------------
diff --git a/src/geodef/plot.py b/src/geodef/plot.py
index f874d5b..d9353da 100644
--- a/src/geodef/plot.py
+++ b/src/geodef/plot.py
@@ -108,18 +108,12 @@ def _stations_to_local_km(
 
     Uses the fault's reference point as the local origin.
     """
-    from geodef import transforms
-
-    alt = np.zeros(dataset.n_stations)
-    e, n, _ = transforms.geod2enu(
-        dataset.lat,
-        dataset.lon,
-        alt,
-        fault._ref_lat,
-        fault._ref_lon,
-        0.0,
+    enu = fault.frame.to_enu(
+        lon=dataset.lon,
+        lat=dataset.lat,
+        alt=np.full(dataset.n_stations, fault.frame.origin_alt),
     )
-    return e * 1e-3, n * 1e-3
+    return enu[:, 0] * 1e-3, enu[:, 1] * 1e-3
 
 
 def _get_patch_vertices_local(fault: Fault) -> list[np.ndarray]:
@@ -132,11 +126,6 @@ def _get_patch_vertices_local(fault: Fault) -> list[np.ndarray]:
         List of arrays, each (n_corners, 2). Rectangular patches have
         4 corners; triangular patches have 3.
     """
-    from geodef import transforms
-
-    ref_lat = fault._ref_lat
-    ref_lon = fault._ref_lon
-
     if fault.engine == "okada":
         assert fault._length is not None and fault._width is not None
         cos_dip = np.cos(np.radians(fault.dip))
@@ -166,16 +155,8 @@ def _get_patch_vertices_local(fault: Fault) -> list[np.ndarray]:
             ]
         )
 
-        # Patch centers in local ENU (meters)
-        alt = np.zeros(fault.n_patches)
-        ce, cn, _ = transforms.geod2enu(
-            fault._lat,
-            fault._lon,
-            alt,
-            ref_lat,
-            ref_lon,
-            0.0,
-        )
+        centers = fault.centers_local
+        ce, cn = centers[:, 0], centers[:, 1]
 
         verts = []
         for i in range(fault.n_patches):
@@ -206,11 +187,6 @@ def _get_patch_vertices_3d(fault: Fault) -> list[np.ndarray]:
     Returns:
         List of arrays, each (n_corners, 3).
     """
-    from geodef import transforms
-
-    ref_lat = fault._ref_lat
-    ref_lon = fault._ref_lon
-
     if fault.engine == "okada":
         assert fault._length is not None and fault._width is not None
         sin_dip = np.sin(np.radians(fault.dip))
@@ -247,15 +223,8 @@ def _get_patch_vertices_3d(fault: Fault) -> list[np.ndarray]:
             ]
         )
 
-        alt = np.zeros(fault.n_patches)
-        ce, cn, _ = transforms.geod2enu(
-            fault._lat,
-            fault._lon,
-            alt,
-            ref_lat,
-            ref_lon,
-            0.0,
-        )
+        centers = fault.centers_local
+        ce, cn = centers[:, 0], centers[:, 1]
 
         verts = []
         for i in range(fault.n_patches):
@@ -354,11 +323,11 @@ def _get_slip_component(
     extract.
 
     Args:
-        slip: Either a single-component vector of length N, or a blocked
+        slip: A single-component vector of length N, or a blocked
             ``[ss_0..ss_N, ds_0..ds_N]`` vector of length 2*N.
         n_patches: Number of fault patches (N).
-        component: One of ``'strike'``, ``'dip'``, ``'magnitude'``. Only
-            used when *slip* has length 2*N.
+        component: One of ``'strike'``, ``'dip'``, ``'magnitude'``,
+            or ``'magnitude'``.
 
     Returns:
         Array of shape (N,).
@@ -366,15 +335,16 @@ def _get_slip_component(
     Raises:
         ValueError: If *component* is invalid or *slip* has wrong length.
     """
-    if slip.shape[0] == n_patches:
-        return slip
-    if slip.shape[0] != 2 * n_patches:
+    slip_array = np.asarray(slip)
+    if slip_array.shape[0] == n_patches:
+        return slip_array
+    if slip_array.shape[0] != 2 * n_patches:
         raise ValueError(
-            f"slip length {slip.shape[0]} does not match "
+            f"slip length {slip_array.shape[0]} does not match "
             f"n_patches = {n_patches} or 2 * n_patches = {2 * n_patches}"
         )
-    ss = slip[:n_patches]
-    ds = slip[n_patches:]
+    ss = slip_array[:n_patches]
+    ds = slip_array[n_patches:]
     if component == "strike":
         return ss
     if component == "dip":
@@ -634,8 +604,8 @@ def slip(
 
     Args:
         fault: Fault geometry (rectangular or triangular).
-        slip_vector: Slip vector, either length N (single component) or
-            length 2*N (blocked ``[ss_0..ss_N, ds_0..ds_N]``).
+        slip_vector: A length-N single component or a length-2*N blocked
+            ``[ss_0..ss_N, ds_0..ds_N]`` vector.
         ax: Axes to plot on. Creates a new figure if ``None``.
         components: Which component to display when *slip_vector* has length
             2*N. One of ``'strike'``, ``'dip'``, or ``'magnitude'``
@@ -721,7 +691,7 @@ def slip_interpolated(
 
     Args:
         fault: Fault geometry (rectangular or triangular).
-        slip_vector: Slip vector, length N or 2*N (see :func:`slip`).
+        slip_vector: Slip vector (see :func:`slip`).
         ax: Axes to plot on. Creates a new figure if ``None``.
         components: Component to display for a 2*N vector: ``'strike'``,
             ``'dip'``, or ``'magnitude'`` (default).
@@ -743,6 +713,8 @@ def slip_interpolated(
             "strike": "Strike-slip (m)",
             "dip": "Dip-slip (m)",
             "magnitude": "Slip magnitude (m)",
+            "rake_parallel": "Rake-parallel slip (m)",
+            "rake_perpendicular": "Rake-perpendicular slip (m)",
         }
         colorbar_label = labels.get(components, "Slip (m)")
 
@@ -1535,9 +1507,8 @@ def map(
             overlaid on the map.
         values: Per-patch scalar array (length *n_patches*) to color
             the patches by. Mutually exclusive with ``slip_vector``.
-        slip_vector: Slip vector, either length *n_patches* or blocked
-            ``[ss | ds]`` length *2 × n_patches*. Decomposed via
-            ``components`` when blocked.
+        slip_vector: Length *n_patches* vector or blocked ``[ss | ds]`` length
+            *2 × n_patches*. Decomposed via ``components`` when blocked.
         components: Which slip component to extract when using
             ``slip_vector``. One of ``'magnitude'``, ``'strike'``,
             ``'dip'`` (default ``'magnitude'``). Ignored for
diff --git a/src/geodef/slip.py b/src/geodef/slip.py
new file mode 100644
index 0000000..29196f9
--- /dev/null
+++ b/src/geodef/slip.py
@@ -0,0 +1,253 @@
+"""Slip-vector packing and basis conversion functions."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+import numpy as np
+import numpy.typing as npt
+
+if TYPE_CHECKING:
+    from geodef.fault import Fault
+
+__all__ = [
+    "from_azimuth",
+    "from_plate",
+    "from_rake",
+    "magnitude",
+    "pack",
+    "plate_rake_from_euler",
+    "rake",
+    "to_plate",
+    "unpack",
+]
+
+
+def _as_1d(values: npt.ArrayLike, name: str) -> np.ndarray:
+    """Return a finite one-dimensional float array."""
+    array = np.atleast_1d(np.asarray(values, dtype=float))
+    if array.ndim != 1:
+        raise ValueError(f"{name} must be one-dimensional, got shape {array.shape}")
+    if not np.all(np.isfinite(array)):
+        raise ValueError(f"{name} must contain only finite values")
+    return array
+
+
+def _matching_components(
+    first: npt.ArrayLike,
+    second: npt.ArrayLike,
+    first_name: str,
+    second_name: str,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Validate two component arrays with the same shape."""
+    first_array = _as_1d(first, first_name)
+    second_array = _as_1d(second, second_name)
+    if first_array.shape != second_array.shape:
+        raise ValueError(
+            f"{first_name} and {second_name} must have the same shape, got "
+            f"{first_array.shape} and {second_array.shape}"
+        )
+    return first_array, second_array
+
+
+def _angle(values: npt.ArrayLike, size: int, name: str) -> np.ndarray:
+    """Broadcast a finite scalar or per-patch angle array."""
+    angle = np.asarray(values, dtype=float)
+    try:
+        result = np.broadcast_to(angle, (size,))
+    except ValueError as exc:
+        raise ValueError(f"{name} must be scalar or have shape ({size},)") from exc
+    if not np.all(np.isfinite(result)):
+        raise ValueError(f"{name} must contain only finite values")
+    return result
+
+
+def pack(strike_slip: npt.ArrayLike, dip_slip: npt.ArrayLike) -> np.ndarray:
+    """Pack physical slip components into GeoDef's blocked vector.
+
+    Args:
+        strike_slip: Strike-slip value per patch, shape ``(N,)``.
+        dip_slip: Dip-slip value per patch, shape ``(N,)``.
+
+    Returns:
+        Blocked ``[strike_slip | dip_slip]`` vector, shape ``(2N,)``.
+
+    Raises:
+        ValueError: If either component is non-finite, not one-dimensional, or
+            has a different shape.
+    """
+    strike_array, dip_array = _matching_components(
+        strike_slip, dip_slip, "strike_slip", "dip_slip"
+    )
+    return np.concatenate([strike_array, dip_array])
+
+
+def unpack(vector: npt.ArrayLike) -> tuple[np.ndarray, np.ndarray]:
+    """Unpack a blocked two-component slip vector.
+
+    Args:
+        vector: Blocked ``[strike_slip | dip_slip]`` vector, shape ``(2N,)``.
+
+    Returns:
+        ``(strike_slip, dip_slip)`` arrays, each shape ``(N,)``.
+
+    Raises:
+        ValueError: If ``vector`` is empty, non-finite, not one-dimensional, or
+            does not contain an even number of entries.
+    """
+    array = _as_1d(vector, "vector")
+    if array.size == 0 or array.size % 2:
+        raise ValueError("vector must contain a non-empty even number of entries")
+    midpoint = array.size // 2
+    return array[:midpoint], array[midpoint:]
+
+
+def from_rake(
+    amplitude: npt.ArrayLike, rake_degrees: npt.ArrayLike
+) -> tuple[np.ndarray, np.ndarray]:
+    """Convert signed amplitudes along rake to strike/dip components.
+
+    Args:
+        amplitude: Signed slip amplitude per patch, shape ``(N,)``.
+        rake_degrees: Local rake in degrees, scalar or shape ``(N,)``.
+
+    Returns:
+        ``(strike_slip, dip_slip)`` arrays, each shape ``(N,)``.
+    """
+    amplitude_array = _as_1d(amplitude, "amplitude")
+    angle = np.deg2rad(_angle(rake_degrees, amplitude_array.size, "rake_degrees"))
+    return amplitude_array * np.cos(angle), amplitude_array * np.sin(angle)
+
+
+def from_azimuth(
+    amplitude: npt.ArrayLike,
+    azimuth_degrees: float,
+    fault_strike_degrees: npt.ArrayLike,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Convert slip along a geographic azimuth to patch-local components.
+
+    Args:
+        amplitude: Signed slip amplitude per patch, shape ``(N,)``.
+        azimuth_degrees: Geographic direction clockwise from North.
+        fault_strike_degrees: Strike of each patch in degrees, shape ``(N,)``.
+
+    Returns:
+        ``(strike_slip, dip_slip)`` arrays, each shape ``(N,)``.
+    """
+    amplitude_array = _as_1d(amplitude, "amplitude")
+    strike = _angle(fault_strike_degrees, amplitude_array.size, "fault_strike_degrees")
+    return from_rake(amplitude_array, float(azimuth_degrees) - strike)
+
+
+def from_plate(
+    parallel: npt.ArrayLike,
+    perpendicular: npt.ArrayLike,
+    plate_rake_degrees: npt.ArrayLike,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Rotate plate-coordinate slip into patch-local strike/dip components.
+
+    Args:
+        parallel: Plate-rake-parallel slip per patch, shape ``(N,)``.
+        perpendicular: Plate-rake-perpendicular slip per patch, shape ``(N,)``.
+        plate_rake_degrees: Plate direction in local rake coordinates, scalar or
+            shape ``(N,)``.
+
+    Returns:
+        ``(strike_slip, dip_slip)`` arrays, each shape ``(N,)``.
+    """
+    parallel_array, perpendicular_array = _matching_components(
+        parallel, perpendicular, "parallel", "perpendicular"
+    )
+    angle = np.deg2rad(
+        _angle(
+            plate_rake_degrees,
+            parallel_array.size,
+            "plate_rake_degrees",
+        )
+    )
+    strike_slip = parallel_array * np.cos(angle) - perpendicular_array * np.sin(angle)
+    dip_slip = parallel_array * np.sin(angle) + perpendicular_array * np.cos(angle)
+    return strike_slip, dip_slip
+
+
+def to_plate(
+    strike_slip: npt.ArrayLike,
+    dip_slip: npt.ArrayLike,
+    plate_rake_degrees: npt.ArrayLike,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Rotate patch-local strike/dip slip into plate coordinates.
+
+    Args:
+        strike_slip: Strike-slip value per patch, shape ``(N,)``.
+        dip_slip: Dip-slip value per patch, shape ``(N,)``.
+        plate_rake_degrees: Plate direction in local rake coordinates, scalar or
+            shape ``(N,)``.
+
+    Returns:
+        ``(parallel, perpendicular)`` arrays, each shape ``(N,)``.
+    """
+    strike_array, dip_array = _matching_components(
+        strike_slip, dip_slip, "strike_slip", "dip_slip"
+    )
+    angle = np.deg2rad(
+        _angle(plate_rake_degrees, strike_array.size, "plate_rake_degrees")
+    )
+    parallel = strike_array * np.cos(angle) + dip_array * np.sin(angle)
+    perpendicular = -strike_array * np.sin(angle) + dip_array * np.cos(angle)
+    return parallel, perpendicular
+
+
+def magnitude(strike_slip: npt.ArrayLike, dip_slip: npt.ArrayLike) -> np.ndarray:
+    """Return unsigned physical slip magnitude per patch.
+
+    Args:
+        strike_slip: Strike-slip value per patch, shape ``(N,)``.
+        dip_slip: Dip-slip value per patch, shape ``(N,)``.
+
+    Returns:
+        Slip magnitude, shape ``(N,)``.
+    """
+    strike_array, dip_array = _matching_components(
+        strike_slip, dip_slip, "strike_slip", "dip_slip"
+    )
+    return np.hypot(strike_array, dip_array)
+
+
+def rake(strike_slip: npt.ArrayLike, dip_slip: npt.ArrayLike) -> np.ndarray:
+    """Return physical local rake in degrees per patch.
+
+    Args:
+        strike_slip: Strike-slip value per patch, shape ``(N,)``.
+        dip_slip: Dip-slip value per patch, shape ``(N,)``.
+
+    Returns:
+        Rake in degrees, shape ``(N,)``.
+    """
+    strike_array, dip_array = _matching_components(
+        strike_slip, dip_slip, "strike_slip", "dip_slip"
+    )
+    return np.degrees(np.arctan2(dip_array, strike_array))
+
+
+def plate_rake_from_euler(fault: Fault, pole: tuple[float, float, float]) -> np.ndarray:
+    """Compute plate direction in each patch's local rake coordinates.
+
+    Args:
+        fault: Fault whose patch centers and strikes define the basis.
+        pole: ``(latitude, longitude, rate)`` in degrees, degrees, deg/Myr.
+
+    Returns:
+        Local plate rake in degrees for each patch, shape ``(N,)``.
+
+    Raises:
+        ValueError: If a patch center lies on the Euler axis and therefore has
+            no defined velocity direction.
+    """
+    from geodef.euler import pole_velocity
+
+    centers = fault.centers_geo
+    east, north = pole_velocity(centers[:, 1], centers[:, 0], pole[0], pole[1], pole[2])
+    if np.any(np.hypot(east, north) == 0.0):
+        raise ValueError("Euler pole produces zero velocity at a patch center")
+    azimuth = np.degrees(np.arctan2(east, north))
+    return np.asarray(azimuth - fault.strike)
diff --git a/tests/test_bayes.py b/tests/test_bayes.py
index 81e8454..fcd6064 100644
--- a/tests/test_bayes.py
+++ b/tests/test_bayes.py
@@ -13,6 +13,7 @@
 from geodef import backend, bayes, gradients
 from geodef.data import GNSS
 from geodef.fault import Fault
+from geodef.geometry import LocalFrame
 from geodef.invert import LinearSystem, _projection_matrix
 
 jax = pytest.importorskip("jax")
@@ -114,6 +115,25 @@ def _posterior(gnss_data, **overrides):
 
 
 class TestConstruction:
+    def test_accepts_parameter_mapping(self, gnss_data):
+        parameters = dict(
+            zip(
+                ["e0", "n0", "depth", "strike", "dip", "length", "width"],
+                _THETA_TRUE,
+                strict=True,
+            )
+        )
+
+        post = _posterior(gnss_data, theta0=parameters)
+
+        assert post.frame == LocalFrame(_REF_LAT, _REF_LON)
+        assert post.geometry(post.x0) == parameters
+        assert post.fault(post.x0).frame == post.frame
+
+    def test_rejects_incompatible_explicit_and_legacy_frames(self, gnss_data):
+        with pytest.raises(ValueError, match="incompatible local frames"):
+            _posterior(gnss_data, frame=LocalFrame(_REF_LAT, _REF_LON + 1.0))
+
     def test_param_names_hierarchical(self, gnss_data):
         post = _posterior(gnss_data)
         assert post.param_names == ["dip", "depth", "log10_sigma", "log10_lambda"]
diff --git a/tests/test_bayes_slip.py b/tests/test_bayes_slip.py
index e21f02a..d9fbec6 100644
--- a/tests/test_bayes_slip.py
+++ b/tests/test_bayes_slip.py
@@ -399,7 +399,7 @@ def _rect_style_collapsed_logpdf(post, log10_sigma, log10_lambda):
     comparison isolates the sigma/lambda power convention from an
     unrelated discrepancy: ``RectPosterior`` assembles G with the
     JAX-autodiff-friendly flat-Cartesian ``rect_greens``, while
-    ``SlipPosterior`` uses the full geodetic ``greens()`` pipeline via
+    ``SlipPosterior`` uses the full geodetic ``matrix()`` pipeline via
     ``LinearSystem`` — the two differ at the ~1e-4 relative level (test
     ``test_bayes.py`` sidesteps the same discrepancy by injecting a
     shared G into a ``LinearSystem`` before comparing to ABIC), enough
diff --git a/tests/test_bayes_tri.py b/tests/test_bayes_tri.py
index 147505a..8a3498d 100644
--- a/tests/test_bayes_tri.py
+++ b/tests/test_bayes_tri.py
@@ -183,6 +183,13 @@ def tri_post_profiled(warp4, gnss_tri):
 
 
 class TestTriWarp:
+    def test_preserves_fault_frame(self, small_mesh_fault):
+        warp = bayes.TriWarp(small_mesh_fault, n_knots=(2, 2))
+
+        assert warp.frame is small_mesh_fault.frame
+        trial = warp.fault(np.zeros(warp.n_knots))
+        assert trial.frame is small_mesh_fault.frame
+
     def test_n_knots_and_shapes(self, small_warp):
         assert small_warp.n_knots == 6
         assert small_warp.knots_uv.shape == (6, 2)
@@ -287,6 +294,9 @@ def test_explicit_knots(self, small_mesh_fault):
 
 
 class TestFrameAnchor:
+    def test_exposes_warp_frame(self, tri_post_hier, warp4):
+        assert tri_post_hier.frame is warp4.frame
+
     def test_matches_linear_system_g_w_at_theta_zero(
         self, tri_post_hier, mesh_fault, gnss_tri
     ):
diff --git a/tests/test_cache.py b/tests/test_cache.py
index a088294..f36e004 100644
--- a/tests/test_cache.py
+++ b/tests/test_cache.py
@@ -9,7 +9,8 @@
 from geodef.cache import compute_hash
 from geodef.data import GNSS, InSAR
 from geodef.fault import Fault
-from geodef.greens import greens
+from geodef.geometry import LocalFrame
+from geodef.greens import matrix as greens
 
 # ====================================================================
 # Group 1: compute_hash (pure function)
@@ -296,7 +297,7 @@ def test_info_after_writes(self, tmp_path: Path) -> None:
 
 
 # ====================================================================
-# Group 5: Integration with greens()
+# Group 5: Integration with matrix()
 # ====================================================================
 
 
@@ -352,7 +353,7 @@ def insar_data() -> InSAR:
 
 
 class TestGreensIntegration:
-    """Tests for caching integration with greens()."""
+    """Tests for caching integration with matrix()."""
 
     def test_greens_caches_result(self, fault_small: Fault, gnss_data: GNSS) -> None:
         """First call creates cache file, second call uses it."""
@@ -391,6 +392,25 @@ def test_greens_invalidates_on_fault_change(self, gnss_data: GNSS) -> None:
         greens(fault_b, gnss_data)
         assert cache.info()["n_files"] == 2
 
+    def test_tri_greens_invalidates_on_frame_change(self, gnss_data: GNSS) -> None:
+        """A triangular vertex frame is a numerical input to Green's assembly."""
+        vertices = np.array(
+            [
+                [
+                    [0.0, 0.0, -10_000.0],
+                    [10_000.0, 0.0, -10_000.0],
+                    [0.0, 10_000.0, -15_000.0],
+                ]
+            ]
+        )
+        fault_a = Fault.from_triangles(vertices, frame=LocalFrame(0.0, 100.0))
+        fault_b = Fault.from_triangles(vertices, frame=LocalFrame(0.0, 101.0))
+
+        greens(fault_a, gnss_data)
+        greens(fault_b, gnss_data)
+
+        assert cache.info()["n_files"] == 2
+
     def test_greens_invalidates_on_data_change(self, fault_small: Fault) -> None:
         """Different observation locations produce a different cache entry."""
         lat_a = np.array([0.2, -0.2])
diff --git a/tests/test_conventions.py b/tests/test_conventions.py
index 4cc7c9d..395d30c 100644
--- a/tests/test_conventions.py
+++ b/tests/test_conventions.py
@@ -24,7 +24,8 @@
 
 import geodef
 from geodef import GNSS, Fault
-from geodef.greens import greens, stack_obs, stack_weights
+from geodef.greens import matrix as greens
+from geodef.greens import stack_obs, stack_weights
 
 LAM = 3.0e-2
 SCALE = 7.5  # arbitrary c != 1
@@ -85,7 +86,7 @@ class TestLinearPaths:
     def test_direct_solve_matches_normal_equations(self, problem) -> None:
         """invert() solves (GtWG + lambda LtL) m = GtWd — lambda, not lambda^2."""
         fault, gnss, L = problem
-        result = geodef.invert(fault, gnss, smoothing=L, smoothing_strength=LAM)
+        result = geodef.invert.solve(fault, gnss, smoothing=L, smoothing_strength=LAM)
         G = greens(fault, gnss)
         W = stack_weights(gnss)
         d = stack_obs(gnss)
@@ -94,7 +95,7 @@ def test_direct_solve_matches_normal_equations(self, problem) -> None:
 
     def test_linear_system_matches_invert(self, problem) -> None:
         fault, gnss, L = problem
-        a = geodef.invert(fault, gnss, smoothing=L, smoothing_strength=LAM)
+        a = geodef.invert.solve(fault, gnss, smoothing=L, smoothing_strength=LAM)
         b = geodef.LinearSystem(fault, gnss, L).invert(smoothing_strength=LAM)
         npt.assert_allclose(a.slip_vector, b.slip_vector, rtol=1e-12)
 
@@ -108,14 +109,14 @@ def test_augmented_system_equivalence(self, problem) -> None:
         G_aug = np.vstack([sqrtW @ G, np.sqrt(LAM) * L])
         d_aug = np.concatenate([sqrtW @ d, np.zeros(L.shape[0])])
         m_aug, *_ = np.linalg.lstsq(G_aug, d_aug, rcond=None)
-        result = geodef.invert(fault, gnss, smoothing=L, smoothing_strength=LAM)
+        result = geodef.invert.solve(fault, gnss, smoothing=L, smoothing_strength=LAM)
         npt.assert_allclose(result.slip_vector, m_aug, rtol=1e-6, atol=1e-10)
 
     def test_invert_scaling_invariance(self, problem) -> None:
         """(L, lam) -> (L/sqrt(c), c lam) must not change the solution."""
         fault, gnss, L = problem
-        a = geodef.invert(fault, gnss, smoothing=L, smoothing_strength=LAM)
-        b = geodef.invert(
+        a = geodef.invert.solve(fault, gnss, smoothing=L, smoothing_strength=LAM)
+        b = geodef.invert.solve(
             fault,
             gnss,
             smoothing=L / np.sqrt(SCALE),
@@ -126,7 +127,7 @@ def test_invert_scaling_invariance(self, problem) -> None:
     def test_model_covariance_convention(self, problem) -> None:
         """C_m = (GtWG + lambda LtL)^-1 with lambda to the first power."""
         fault, gnss, L = problem
-        result = geodef.invert(fault, gnss, smoothing=L, smoothing_strength=LAM)
+        result = geodef.invert.solve(fault, gnss, smoothing=L, smoothing_strength=LAM)
         C = geodef.model_covariance(result, fault, gnss)
         G = greens(fault, gnss)
         W = stack_weights(gnss)
diff --git a/tests/test_docs.py b/tests/test_docs.py
index 0accd80..749bfe5 100644
--- a/tests/test_docs.py
+++ b/tests/test_docs.py
@@ -25,6 +25,7 @@
     "euler",
     "fault",
     "geomap",
+    "geometry",
     "gradients",
     "greens",
     "invert",
diff --git a/tests/test_fault.py b/tests/test_fault.py
index e6d4254..0f610b7 100644
--- a/tests/test_fault.py
+++ b/tests/test_fault.py
@@ -8,6 +8,7 @@
 from pathlib import Path
 
 import numpy as np
+import numpy.testing as npt
 import pytest
 
 from geodef.fault import (
@@ -16,6 +17,8 @@
     magnitude_to_moment,
     moment_to_magnitude,
 )
+from geodef.geometry import LocalFrame
+from geodef.medium import ElasticMedium
 
 # ======================================================================
 # Fixtures
@@ -117,12 +120,73 @@ def test_arrays_are_read_only(self, simple_fault):
         with pytest.raises(ValueError):
             simple_fault._lat[0] = 999.0
 
+    def test_direct_fault_has_explicit_frame(self, single_patch):
+        assert single_patch.frame.projection == "wgs84-enu"
+        assert single_patch.frame.origin_lat == pytest.approx(0.0)
+        assert single_patch.frame.origin_lon == pytest.approx(100.0)
+
 
 # ======================================================================
 # 2. Fault.planar() factory
 # ======================================================================
 
 
+class TestPlanarFrame:
+    """Fault.planar owns its optional local frame without another wrapper."""
+
+    def test_scalar_call_uses_centered_frame_by_default(self):
+        fault = Fault.planar(
+            lat=-2.0,
+            lon=100.0,
+            depth=15_000.0,
+            strike=315.0,
+            dip=25.0,
+            length=80_000.0,
+            width=40_000.0,
+        )
+
+        assert fault.frame.origin_lat == -2.0
+        assert fault.frame.origin_lon == 100.0
+        npt.assert_allclose(np.mean(fault.centers_local, axis=0)[:2], 0.0, atol=1.0)
+
+    def test_explicit_frame_survives_medium_copy(self):
+        frame = LocalFrame(0.0, 100.0)
+        fault = Fault.planar(
+            lat=0.01,
+            lon=100.02,
+            depth=10_000.0,
+            strike=0.0,
+            dip=30.0,
+            length=20_000.0,
+            width=10_000.0,
+            frame=frame,
+        )
+
+        changed = fault.with_medium(ElasticMedium(shear_modulus=40e9))
+
+        assert changed.frame is frame
+        npt.assert_allclose(changed.centers_geo, fault.centers_geo)
+
+    def test_to_frame_preserves_geographic_geometry(self):
+        fault = Fault.planar(
+            lat=0.0,
+            lon=100.0,
+            depth=10_000.0,
+            strike=0.0,
+            dip=30.0,
+            length=20_000.0,
+            width=10_000.0,
+            n_length=2,
+            n_width=2,
+        )
+        target = LocalFrame(0.1, 100.2)
+
+        transformed = fault.to_frame(target)
+
+        assert transformed.frame is target
+        npt.assert_allclose(transformed.centers_geo, fault.centers_geo)
+
+
 class TestPlanar:
     """Test Fault.planar() factory classmethod."""
 
@@ -395,6 +459,37 @@ def test_patch_index_no_grid_raises(self):
         with pytest.raises(ValueError, match="structured grid"):
             fault.patch_index(0, 0)
 
+    def test_reshape_and_flatten_patches_roundtrip(self, simple_fault):
+        values = np.arange(simple_fault.n_patches)
+
+        grid = simple_fault.reshape_patches(values)
+
+        assert grid.shape == (5, 10)
+        np.testing.assert_array_equal(grid[0], np.arange(10))
+        np.testing.assert_array_equal(simple_fault.flatten_patches(grid), values)
+
+    def test_reshape_and_flatten_preserve_trailing_dimensions(self, simple_fault):
+        values = np.arange(simple_fault.n_patches * 2).reshape(-1, 2)
+
+        grid = simple_fault.reshape_patches(values)
+
+        assert grid.shape == (5, 10, 2)
+        np.testing.assert_array_equal(simple_fault.flatten_patches(grid), values)
+
+    def test_reshape_patches_rejects_unstructured_fault(self):
+        fault = Fault(
+            np.array([0.0]),
+            np.array([100.0]),
+            np.array([10e3]),
+            np.array([0.0]),
+            np.array([90.0]),
+            np.array([10e3]),
+            np.array([10e3]),
+            grid_shape=None,
+        )
+        with pytest.raises(ValueError, match="structured grid"):
+            fault.reshape_patches([1.0])
+
 
 # ======================================================================
 # 8. File I/O
diff --git a/tests/test_geometry.py b/tests/test_geometry.py
new file mode 100644
index 0000000..7b77c6b
--- /dev/null
+++ b/tests/test_geometry.py
@@ -0,0 +1,177 @@
+"""Tests for coordinate frames and geometry-array functions."""
+
+from dataclasses import FrozenInstanceError
+
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+from geodef.geometry import (
+    LocalFrame,
+    as_planar_vector,
+    planar_parameter_dict,
+    triangle_strike_dip,
+    vertices_from_nodes,
+)
+
+
+class TestLocalFrame:
+    """LocalFrame validation, conversion, and compatibility."""
+
+    def test_records_projection_and_is_immutable(self) -> None:
+        frame = LocalFrame(1.0, 100.0, origin_alt=12.0)
+
+        assert frame.projection == "wgs84-enu"
+        with pytest.raises(FrozenInstanceError):
+            frame.origin_lat = 2.0  # type: ignore[misc]
+
+    @pytest.mark.parametrize(
+        ("kwargs", "match"),
+        [
+            ({"origin_lat": 91.0, "origin_lon": 0.0}, "origin_lat"),
+            ({"origin_lat": 0.0, "origin_lon": np.nan}, "origin_lon"),
+            (
+                {
+                    "origin_lat": 0.0,
+                    "origin_lon": 0.0,
+                    "projection": "utm",
+                },
+                "projection",
+            ),
+        ],
+    )
+    def test_invalid_definition_raises(
+        self, kwargs: dict[str, float | str], match: str
+    ) -> None:
+        with pytest.raises(ValueError, match=match):
+            LocalFrame(**kwargs)  # type: ignore[arg-type]
+
+    def test_geographic_enu_round_trip(self) -> None:
+        frame = LocalFrame(1.0, 100.0, origin_alt=25.0)
+        lon = np.array([100.0, 100.02])
+        lat = np.array([1.0, 1.01])
+        alt = np.array([25.0, -1500.0])
+
+        enu = frame.to_enu(lon=lon, lat=lat, alt=alt)
+        geographic = frame.to_geographic(east=enu[:, 0], north=enu[:, 1], up=enu[:, 2])
+
+        assert enu.shape == (2, 3)
+        npt.assert_allclose(geographic[:, 0], lon, atol=1e-10)
+        npt.assert_allclose(geographic[:, 1], lat, atol=1e-10)
+        npt.assert_allclose(geographic[:, 2], alt, atol=1e-6)
+
+    def test_explicit_transform_between_frames(self) -> None:
+        source = LocalFrame(0.0, 100.0)
+        target = LocalFrame(0.1, 100.2)
+        coordinates = np.array([[0.0, 0.0, -1000.0], [500.0, 200.0, 0.0]])
+
+        transformed = source.transform_enu(coordinates, target=target)
+        geographic = source.to_geographic(
+            east=coordinates[:, 0],
+            north=coordinates[:, 1],
+            up=coordinates[:, 2],
+        )
+        expected = target.to_enu(
+            lon=geographic[:, 0],
+            lat=geographic[:, 1],
+            alt=geographic[:, 2],
+        )
+
+        npt.assert_allclose(transformed, expected)
+
+    def test_incompatible_frames_are_rejected(self) -> None:
+        frame = LocalFrame(0.0, 100.0)
+        other = LocalFrame(0.0, 101.0)
+
+        assert frame.is_compatible(frame)
+        assert not frame.is_compatible(other)
+        with pytest.raises(ValueError, match="incompatible local frames"):
+            frame.require_compatible(other)
+
+
+class TestPlanarFunctions:
+    """Planar parameter mappings and expert vectors round-trip."""
+
+    def test_mapping_to_vector_and_back(self) -> None:
+        parameters = {
+            "depth": 15_000.0,
+            "e0": 1200.0,
+            "n0": -500.0,
+            "strike": 315.0,
+            "dip": 25.0,
+            "length": 80_000.0,
+            "width": 40_000.0,
+        }
+
+        vector = as_planar_vector(parameters)
+
+        npt.assert_array_equal(
+            vector,
+            [1200.0, -500.0, 15_000.0, 315.0, 25.0, 80_000.0, 40_000.0],
+        )
+        assert planar_parameter_dict(vector) == parameters
+
+    @pytest.mark.parametrize(
+        ("field", "value"),
+        [
+            ("depth", -1.0),
+            ("strike", 360.0),
+            ("dip", 91.0),
+            ("length", 0.0),
+            ("width", np.inf),
+        ],
+    )
+    def test_invalid_parameter_raises(self, field: str, value: float) -> None:
+        parameters = {
+            "e0": 0.0,
+            "n0": 0.0,
+            "depth": 10_000.0,
+            "strike": 0.0,
+            "dip": 30.0,
+            "length": 20_000.0,
+            "width": 10_000.0,
+        }
+        parameters[field] = value
+
+        with pytest.raises(ValueError, match=field):
+            as_planar_vector(parameters)
+
+
+class TestTriangleFunctions:
+    """Triangle expansion and orientation functions."""
+
+    def test_orientation(self) -> None:
+        vertices = np.array(
+            [
+                [[0.0, 0.0, -1000.0], [1000.0, 0.0, -1000.0], [0.0, 0.0, -2000.0]],
+                [
+                    [1000.0, 0.0, -1000.0],
+                    [1000.0, 0.0, -2000.0],
+                    [0.0, 0.0, -2000.0],
+                ],
+            ]
+        )
+
+        strike, dip = triangle_strike_dip(vertices)
+
+        npt.assert_allclose(dip, 90.0)
+        assert np.all((strike >= 0.0) & (strike < 360.0))
+
+    def test_vertices_from_nodes_preserves_connectivity_order(self) -> None:
+        nodes = np.array(
+            [
+                [0.0, 0.0, -1000.0],
+                [1000.0, 0.0, -1000.0],
+                [0.0, 1000.0, -2000.0],
+                [1000.0, 1000.0, -2000.0],
+            ]
+        )
+        triangles = np.array([[1, 3, 2], [0, 1, 2]])
+
+        vertices = vertices_from_nodes(nodes, triangles)
+
+        npt.assert_array_equal(vertices[0], nodes[triangles[0]])
+
+    def test_vertices_from_nodes_rejects_invalid_connectivity(self) -> None:
+        with pytest.raises(ValueError, match="triangles"):
+            vertices_from_nodes(np.zeros((3, 3)), np.array([[0, 1, 3]]))
diff --git a/tests/test_geometry_search.py b/tests/test_geometry_search.py
index ffe7ee7..2c570c1 100644
--- a/tests/test_geometry_search.py
+++ b/tests/test_geometry_search.py
@@ -11,6 +11,7 @@
 from geodef import backend
 from geodef.data import GNSS
 from geodef.fault import Fault
+from geodef.geometry import LocalFrame
 from geodef.invert import geometry_search
 
 jax = pytest.importorskip("jax")
@@ -108,6 +109,40 @@ def _theta_start(**overrides):
 
 
 class TestGeometrySearch:
+    def test_accepts_mapping_and_returns_fault(self, gnss_data):
+        frame = LocalFrame(_REF_LAT, _REF_LON)
+        parameters = dict(
+            zip(
+                ["e0", "n0", "depth", "strike", "dip", "length", "width"],
+                _theta_start(dip=30.0),
+                strict=True,
+            )
+        )
+
+        result = geometry_search(
+            parameters,
+            gnss_data,
+            frame=frame,
+            free=["dip"],
+            bounds={"dip": (5.0, 45.0)},
+            **_THETA0_KWARGS,
+        )
+
+        assert isinstance(result.fault, Fault)
+        assert result.fault.frame is frame
+        assert result.frame is frame
+        assert abs(np.mean(result.fault.dip) - _TRUE["dip"]) < 0.5
+
+    def test_rejects_mapping_with_missing_parameter(self, gnss_data):
+        with pytest.raises(ValueError, match="missing keys"):
+            geometry_search(
+                {"depth": 20_000.0},
+                gnss_data,
+                frame=LocalFrame(_REF_LAT, _REF_LON),
+                free=["dip"],
+                **_THETA0_KWARGS,
+            )
+
     def test_recovers_dip(self, gnss_data):
         result = geometry_search(
             _theta_start(dip=30.0),
diff --git a/tests/test_greens.py b/tests/test_greens.py
index ac84896..06c3a5a 100644
--- a/tests/test_greens.py
+++ b/tests/test_greens.py
@@ -10,8 +10,35 @@
     build_laplacian_2d,
     build_laplacian_2d_simple,
     build_laplacian_knn,
+    laplacian,
+    select_slip_columns,
 )
 
+
+class TestSlipColumnSelection:
+    """Green's matrices rotate into declared model coordinates."""
+
+    def test_plate_rake_keeps_parallel_and_perpendicular_columns(self):
+        n = 2
+        strike = np.array([[1.0, 2.0], [3.0, 4.0]])
+        dip = np.array([[5.0, 6.0], [7.0, 8.0]])
+        full = np.hstack([strike, dip])
+
+        selected = select_slip_columns(
+            full, n, "plate", plate_rake=np.array([0.0, 90.0])
+        )
+
+        expected_parallel = np.array([[1.0, 6.0], [3.0, 8.0]])
+        expected_perpendicular = np.array([[5.0, -2.0], [7.0, -4.0]])
+        np.testing.assert_allclose(
+            selected, np.hstack([expected_parallel, expected_perpendicular])
+        )
+
+    def test_plate_rake_is_required(self):
+        with pytest.raises(ValueError, match="plate_rake"):
+            select_slip_columns(np.ones((3, 4)), 2, "plate")
+
+
 # ---------------------------------------------------------------------------
 # Laplacian (forward/backward difference boundaries)
 # ---------------------------------------------------------------------------
@@ -20,6 +47,23 @@
 class TestLaplacian2D:
     """Tests for build_laplacian_2d."""
 
+    def test_fault_laplacian_function(self):
+        from geodef.fault import Fault
+
+        fault = Fault.planar(
+            lat=0.0,
+            lon=100.0,
+            depth=10_000.0,
+            strike=0.0,
+            dip=30.0,
+            length=20_000.0,
+            width=10_000.0,
+            n_length=3,
+            n_width=3,
+        )
+
+        np.testing.assert_array_equal(laplacian(fault), fault.laplacian)
+
     def test_matrix_shape(self):
         nL, nW = 5, 4
         L = build_laplacian_2d(nL, nW)
diff --git a/tests/test_greens_integration.py b/tests/test_greens_integration.py
index 644bb87..46fbd31 100644
--- a/tests/test_greens_integration.py
+++ b/tests/test_greens_integration.py
@@ -1,8 +1,8 @@
-"""Tests for geodef.greens() polymorphic Green's matrix assembly (Phase 3.3).
+"""Tests for polymorphic Green's matrix assembly (Phase 3.3).
 
-Covers: greens() with single and joint datasets, all data types (GNSS,
+Covers: matrix() with single and joint datasets, all data types (GNSS,
 InSAR, Vertical), consistency with fault.displacement(), stack_obs(),
-stack_weights(), _project_greens(), and tri engine support.
+stack_weights(), project(), and tri engine support.
 """
 
 import numpy as np
@@ -12,8 +12,10 @@
 from geodef.data import GNSS, InSAR, Vertical
 from geodef.fault import Fault
 from geodef.greens import (
-    _project_greens,
-    greens,
+    matrix as greens,
+)
+from geodef.greens import (
+    project,
     select_slip_columns,
     stack_obs,
     stack_weights,
@@ -123,7 +125,7 @@ def vertical_4pt(obs_points):
 
 
 # ======================================================================
-# 1. greens() shape tests — single dataset
+# 1. matrix() shape tests — single dataset
 # ======================================================================
 
 
@@ -148,7 +150,7 @@ def test_vertical_shape(self, fault_4x3, vertical_4pt):
 
 
 # ======================================================================
-# 2. greens() shape tests — joint datasets
+# 2. matrix() shape tests — joint datasets
 # ======================================================================
 
 
@@ -180,7 +182,7 @@ def test_joint_equals_individual_vstack(
 
 
 class TestConsistencyWithDisplacement:
-    """Verify that greens() @ slip == fault.displacement() results."""
+    """Verify that matrix() @ slip == fault.displacement() results."""
 
     def test_gnss_forward_model(self, single_patch, obs_points):
         lat, lon = obs_points
@@ -360,21 +362,21 @@ def test_stack_weights_positive_definite(self, gnss_4station, insar_4pixel):
 
 
 # ======================================================================
-# 7. _project_greens internal helper
+# 7. project() helper
 # ======================================================================
 
 
 class TestProjectGreens:
-    """Tests for the _project_greens helper function."""
+    """Tests for the project helper function."""
 
     def test_insar_projection_reduces_rows(self, fault_4x3, insar_4pixel):
         G_raw = fault_4x3.greens_matrix(insar_4pixel.lat, insar_4pixel.lon)
-        G_proj = _project_greens(insar_4pixel, G_raw)
+        G_proj = project(insar_4pixel, G_raw)
         assert G_proj.shape == (4, 24)  # from (12, 24) to (4, 24)
 
     def test_gnss_3comp_preserves_rows(self, fault_4x3, gnss_4station):
         G_raw = fault_4x3.greens_matrix(gnss_4station.lat, gnss_4station.lon)
-        G_proj = _project_greens(gnss_4station, G_raw)
+        G_proj = project(gnss_4station, G_raw)
         assert G_proj.shape == (12, 24)  # 3*4 stays 3*4
 
     def test_vertical_extracts_uz(self, single_patch, obs_points):
@@ -382,7 +384,7 @@ def test_vertical_extracts_uz(self, single_patch, obs_points):
         n = len(lat)
         vert = Vertical(lon=lon, lat=lat, displacement=np.zeros(n), sigma=np.ones(n))
         G_raw = single_patch.greens_matrix(lat, lon)
-        G_proj = _project_greens(vert, G_raw)
+        G_proj = project(vert, G_raw)
 
         # Vertical projection should extract rows 2, 5, 8, 11 from G_raw
         G_uz = G_raw[2::3, :]
@@ -395,10 +397,10 @@ def test_vertical_extracts_uz(self, single_patch, obs_points):
 
 
 class TestTopLevelAPI:
-    """Verify greens() is accessible from geodef namespace."""
+    """Verify matrix() is accessible from the greens namespace."""
 
     def test_greens_accessible(self):
-        assert hasattr(geodef.greens, "greens")
+        assert hasattr(geodef.greens, "matrix")
 
     def test_stack_obs_accessible(self):
         assert hasattr(geodef, "stack_obs")
@@ -469,7 +471,7 @@ def test_tri_greens_nonzero(self, tri_fault):
         assert np.any(G != 0.0)
 
     def test_tri_with_greens_function(self, tri_fault):
-        """greens() works with tri engine faults."""
+        """matrix() works with tri engine faults."""
         obs_lat = np.array([0.1, -0.1])
         obs_lon = np.array([100.1, 99.9])
         n = len(obs_lat)
@@ -501,12 +503,12 @@ def test_tri_invalid_kind_raises(self, tri_fault):
 
 
 # ======================================================================
-# 6. greens() component selection
+# 6. matrix() component selection
 # ======================================================================
 
 
 class TestGreensComponents:
-    """greens(components=...) reduces columns like the inversion path."""
+    """matrix(components=...) reduces columns like the inversion path."""
 
     def test_both_is_default(self, fault_4x3, gnss_4station):
         G_default = greens(fault_4x3, gnss_4station)
diff --git a/tests/test_invert.py b/tests/test_invert.py
index dbdbe9c..92ca1f1 100644
--- a/tests/test_invert.py
+++ b/tests/test_invert.py
@@ -11,16 +11,20 @@
 
 from geodef.data import GNSS, InSAR, Vertical
 from geodef.fault import Fault
-from geodef.greens import greens, stack_obs, stack_weights
+from geodef.greens import matrix as greens
+from geodef.greens import stack_obs, stack_weights
 from geodef.invert import (
     DatasetDiagnostics,
     InversionResult,
     dataset_diagnostics,
-    invert,
     model_covariance,
     model_resolution,
     model_uncertainty,
 )
+from geodef.invert import (
+    solve as invert,
+)
+from geodef.slip import from_plate, pack
 
 # ======================================================================
 # Fixtures
@@ -149,6 +153,15 @@ def test_slip_vector_is_blocked(self, fault_4x3, obs_points):
         np.testing.assert_array_equal(result.slip[:, 0], result.slip_vector[:n])
         np.testing.assert_array_equal(result.slip[:, 1], result.slip_vector[n:])
 
+    def test_result_exposes_named_physical_components(self, fault_4x3, obs_points):
+        gnss = _make_gnss(fault_4x3, obs_points, np.ones(12), np.zeros(12))
+
+        result = invert(fault_4x3, gnss)
+
+        np.testing.assert_allclose(result.strike_slip, result.slip_vector[:12])
+        np.testing.assert_allclose(result.dip_slip, result.slip_vector[12:])
+        np.testing.assert_allclose(result.slip_magnitude, np.abs(result.strike_slip))
+
     def test_residuals_shape(self, fault_4x3, obs_points):
         slip_ss = np.ones(12)
         slip_ds = np.zeros(12)
@@ -1733,6 +1746,25 @@ def test_smoothing_laplacian(self, fault_4x3, obs_points):
         )
         assert result.slip.shape == (12, 1)
 
+    def test_result_converts_one_component_rake(self, fault_4x3, obs_points):
+        rake = 30.0
+        radians = np.deg2rad(rake)
+        gnss = _make_gnss(
+            fault_4x3,
+            obs_points,
+            np.cos(radians) * np.ones(12),
+            np.sin(radians) * np.ones(12),
+        )
+
+        result = invert(fault_4x3, gnss, components="rake", rake=rake)
+
+        np.testing.assert_allclose(
+            result.strike_slip, result.slip_vector * np.cos(radians)
+        )
+        np.testing.assert_allclose(
+            result.dip_slip, result.slip_vector * np.sin(radians)
+        )
+
     def test_save_load_roundtrip(self, fault_4x3, obs_points, tmp_path):
         rake = 30.0
         r = np.deg2rad(rake)
@@ -1770,6 +1802,86 @@ def test_save_table_dip_column_name(self, fault_4x3, obs_points, tmp_path):
         assert "slip_strike_m" not in text
 
 
+# ======================================================================
+# Plate-rake coordinate inversion
+# ======================================================================
+
+
+class TestComponentsPlate:
+    """Plate coordinates smooth and constrain the large-scale basis."""
+
+    def test_recovers_plate_parallel_and_perpendicular(self, fault_4x3, obs_points):
+        plate_rake = np.linspace(20.0, 40.0, fault_4x3.n_patches)
+        truth_strike, truth_dip = from_plate(np.ones(12), np.full(12, 0.2), plate_rake)
+        gnss = _make_gnss(fault_4x3, obs_points, truth_strike, truth_dip)
+
+        result = invert(
+            fault_4x3,
+            gnss,
+            components="plate",
+            plate_rake=plate_rake,
+        )
+
+        assert result.components == "plate"
+        np.testing.assert_allclose(result.plate_rake, plate_rake)
+        np.testing.assert_allclose(result.rake_parallel, 1.0, atol=0.1)
+        np.testing.assert_allclose(result.rake_perpendicular, 0.2, atol=0.1)
+
+    def test_laplacian_and_bounds_use_plate_coordinates(self, fault_4x3, obs_points):
+        plate_rake = np.linspace(20.0, 40.0, fault_4x3.n_patches)
+        truth_strike, truth_dip = from_plate(np.ones(12), np.zeros(12), plate_rake)
+        gnss = _make_gnss(fault_4x3, obs_points, truth_strike, truth_dip)
+
+        result = invert(
+            fault_4x3,
+            gnss,
+            components="plate",
+            plate_rake=plate_rake,
+            smoothing="laplacian",
+            smoothing_strength=1.0,
+            bounds=(np.array([0.0, -0.01]), np.array([2.0, 0.01])),
+        )
+
+        assert np.all(result.rake_parallel >= 0.0)
+        assert np.all(np.abs(result.rake_perpendicular) <= 0.010001)
+
+    def test_plate_vector_is_accepted_as_smoothing_target(self, fault_4x3, obs_points):
+        plate_rake = np.full(12, 30.0)
+        target = pack(np.ones(12), np.zeros(12))
+        truth_strike, truth_dip = from_plate(np.ones(12), np.zeros(12), plate_rake)
+        gnss = _make_gnss(fault_4x3, obs_points, truth_strike, truth_dip)
+
+        result = invert(
+            fault_4x3,
+            gnss,
+            components="plate",
+            plate_rake=plate_rake,
+            smoothing="damping",
+            smoothing_strength=1.0,
+            smoothing_target=target,
+        )
+
+        assert result.components == "plate"
+
+    def test_missing_plate_rake_raises(self, fault_4x3, obs_points):
+        gnss = _make_gnss(fault_4x3, obs_points, np.ones(12), np.zeros(12))
+        with pytest.raises(ValueError, match="plate_rake"):
+            invert(fault_4x3, gnss, components="plate")
+
+    def test_save_load_preserves_plate_basis(self, fault_4x3, obs_points, tmp_path):
+        plate_rake = np.linspace(20.0, 40.0, fault_4x3.n_patches)
+        truth_strike, truth_dip = from_plate(np.ones(12), np.zeros(12), plate_rake)
+        gnss = _make_gnss(fault_4x3, obs_points, truth_strike, truth_dip)
+        result = invert(fault_4x3, gnss, components="plate", plate_rake=plate_rake)
+
+        path = tmp_path / "plate_result.npz"
+        result.save(path)
+        loaded = InversionResult.load(path)
+
+        np.testing.assert_allclose(loaded.slip_vector, result.slip_vector)
+        np.testing.assert_allclose(loaded.plate_rake, plate_rake)
+
+
 # ======================================================================
 # Fixed slip-azimuth inversion
 # ======================================================================
diff --git a/tests/test_linear_system.py b/tests/test_linear_system.py
index cfdf855..60a2f9e 100644
--- a/tests/test_linear_system.py
+++ b/tests/test_linear_system.py
@@ -20,12 +20,14 @@
     LinearSystem,
     abic_curve,
     dataset_diagnostics,
-    invert,
     lcurve,
     model_covariance,
     model_resolution,
     model_uncertainty,
 )
+from geodef.invert import (
+    solve as invert,
+)
 
 # ======================================================================
 # Fixtures (shared with test_invert.py pattern)
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 503c885..55d32bd 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -6,6 +6,8 @@
 import numpy.testing as npt
 import pytest
 
+from geodef.fault import Fault
+from geodef.geometry import LocalFrame
 from geodef.mesh import Mesh, _compute_strike_dip
 
 requires_meshpy = pytest.mark.skipif(
@@ -73,6 +75,34 @@ def test_basic_construction(self, simple_mesh):
         assert simple_mesh.n_nodes == 4
         assert simple_mesh.n_triangles == 2
 
+    def test_records_inferred_frame(self, simple_mesh):
+        assert simple_mesh.frame.projection == "wgs84-enu"
+        assert simple_mesh.frame.origin_lat == pytest.approx(np.mean(simple_mesh.lat))
+        assert simple_mesh.frame.origin_lon == pytest.approx(np.mean(simple_mesh.lon))
+
+    def test_accepts_explicit_frame(self):
+        frame = LocalFrame(0.0, 100.0)
+        mesh = Mesh(
+            lon=np.array([100.0, 100.01, 100.0]),
+            lat=np.array([0.0, 0.0, 0.01]),
+            depth=np.array([0.0, 1000.0, 2000.0]),
+            triangles=np.array([[0, 1, 2]]),
+            frame=frame,
+        )
+
+        assert mesh.frame is frame
+        npt.assert_allclose(mesh.vertices_enu(), mesh.vertices_enu(frame=frame))
+
+    def test_to_frame_preserves_geographic_nodes(self, simple_mesh):
+        frame = LocalFrame(1.0, 101.0)
+
+        transformed = simple_mesh.to_frame(frame)
+
+        assert transformed.frame is frame
+        npt.assert_array_equal(transformed.lon, simple_mesh.lon)
+        npt.assert_array_equal(transformed.lat, simple_mesh.lat)
+        npt.assert_array_equal(transformed.depth, simple_mesh.depth)
+
     def test_arrays_stored(self, simple_mesh):
         assert simple_mesh.lon.shape == (4,)
         assert simple_mesh.lat.shape == (4,)
@@ -113,6 +143,14 @@ def test_shape(self, simple_mesh):
         verts = simple_mesh.vertices_enu(ref_lat=0.0, ref_lon=100.0)
         assert verts.shape == (2, 3, 3)
 
+    def test_rejects_ambiguous_frame_arguments(self, simple_mesh):
+        with pytest.raises(ValueError, match="either frame or ref_lat"):
+            simple_mesh.vertices_enu(
+                ref_lat=0.0,
+                ref_lon=100.0,
+                frame=LocalFrame(0.0, 100.0),
+            )
+
     def test_z_convention(self, simple_mesh):
         """Depth positive down → z negative in ENU (up-positive)."""
         verts = simple_mesh.vertices_enu(ref_lat=0.0, ref_lon=100.0)
@@ -311,6 +349,53 @@ def test_basic_construction(self):
         assert fault.n_patches == 2
         assert fault.engine == "tri"
 
+    def test_accepts_explicit_tri_frame(self):
+        frame = LocalFrame(0.0, 100.0)
+        vertices = np.array(
+            [
+                [
+                    [0.0, 0.0, -1000.0],
+                    [1000.0, 0.0, -1000.0],
+                    [0.0, 0.0, -2000.0],
+                ]
+            ]
+        )
+
+        fault = Fault.from_triangles(vertices, frame=frame)
+
+        assert fault.frame is frame
+        npt.assert_allclose(fault.centers_local, np.mean(vertices, axis=1))
+
+    def test_fault_to_frame_transforms_tri_vertices(self):
+        vertices = np.array(
+            [
+                [
+                    [0.0, 0.0, -1000.0],
+                    [1000.0, 0.0, -1000.0],
+                    [0.0, 0.0, -2000.0],
+                ]
+            ]
+        )
+        fault = Fault.from_triangles(vertices, frame=LocalFrame(0.0, 100.0))
+
+        transformed = fault.to_frame(LocalFrame(0.1, 100.2))
+
+        assert transformed.frame != fault.frame
+        assert transformed.vertices is not None
+        original_geo = fault.frame.to_geographic(
+            east=vertices[..., 0], north=vertices[..., 1], up=vertices[..., 2]
+        )
+        transformed_geo = transformed.frame.to_geographic(
+            east=transformed.vertices[..., 0],
+            north=transformed.vertices[..., 1],
+            up=transformed.vertices[..., 2],
+        )
+        npt.assert_allclose(
+            transformed_geo,
+            original_geo,
+            atol=1e-6,
+        )
+
     def test_strike_dip_derived(self):
         from geodef.fault import Fault
 
@@ -428,6 +513,7 @@ def test_basic(self, simple_mesh):
         fault = Fault.from_mesh(simple_mesh)
         assert fault.n_patches == simple_mesh.n_triangles
         assert fault.engine == "tri"
+        assert fault.frame is simple_mesh.frame
 
     def test_vertices_shape(self, simple_mesh):
         from geodef.fault import Fault
@@ -1393,7 +1479,7 @@ def test_from_trace_to_fault_to_greens(self):
             su=np.array([0.001]),
         )
 
-        G = greens_mod.greens(fault, gnss)
+        G = greens_mod.matrix(fault, gnss)
         assert G.shape[0] == 3  # 3 components
         assert G.shape[1] == fault.n_patches * 2  # ss + ds
 
diff --git a/tests/test_package.py b/tests/test_package.py
index e362dd9..da5b3c8 100644
--- a/tests/test_package.py
+++ b/tests/test_package.py
@@ -4,6 +4,8 @@
 the okada dispatcher works correctly, and the top-level API is usable.
 """
 
+from types import ModuleType
+
 import numpy as np
 import pytest
 
@@ -124,6 +126,16 @@ def test_import_data_module(self):
         assert hasattr(data, "Vertical")
         assert hasattr(data, "DataSet")
 
+    def test_function_oriented_module_namespaces(self):
+        import geodef
+
+        assert isinstance(geodef.invert, ModuleType)
+        assert geodef.solve is geodef.invert.solve
+        assert callable(geodef.greens.matrix)
+        assert callable(geodef.greens.project)
+        assert callable(geodef.greens.laplacian)
+        assert not hasattr(geodef.invert, "invert")
+
 
 # ---------------------------------------------------------------------------
 # 3. Okada dispatcher tests
diff --git a/tests/test_plot.py b/tests/test_plot.py
index d1fb3c3..257c554 100644
--- a/tests/test_plot.py
+++ b/tests/test_plot.py
@@ -310,6 +310,16 @@ def test_rect_fault(self, rect_fault, slip_magnitude):
         ax = geodef.plot.slip(rect_fault, slip_magnitude)
         assert len(ax.collections) >= 1
 
+    def test_blocked_slip_vector(self, rect_fault):
+        vector = np.concatenate(
+            [
+                np.ones(rect_fault.n_patches),
+                np.full(rect_fault.n_patches, 0.5),
+            ]
+        )
+        ax = geodef.plot.slip(rect_fault, vector, components="dip")
+        assert isinstance(ax, plt.Axes)
+
     def test_tri_fault(self, tri_fault, slip_tri):
         ax = geodef.plot.slip(tri_fault, slip_tri)
         assert isinstance(ax, plt.Axes)
diff --git a/tests/test_slip.py b/tests/test_slip.py
new file mode 100644
index 0000000..a136545
--- /dev/null
+++ b/tests/test_slip.py
@@ -0,0 +1,95 @@
+"""Tests for function-oriented slip conversions."""
+
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+from geodef.fault import Fault
+from geodef.slip import (
+    from_azimuth,
+    from_plate,
+    from_rake,
+    magnitude,
+    pack,
+    plate_rake_from_euler,
+    rake,
+    to_plate,
+    unpack,
+)
+
+
+def test_pack_and_unpack_blocked_components() -> None:
+    vector = pack(strike_slip=[1.0, 2.0], dip_slip=[3.0, 4.0])
+
+    strike_slip, dip_slip = unpack(vector)
+
+    npt.assert_allclose(vector, [1.0, 2.0, 3.0, 4.0])
+    npt.assert_allclose(strike_slip, [1.0, 2.0])
+    npt.assert_allclose(dip_slip, [3.0, 4.0])
+
+
+def test_pack_rejects_mismatched_components() -> None:
+    with pytest.raises(ValueError, match="same shape"):
+        pack(strike_slip=[1.0, 2.0], dip_slip=[3.0])
+
+
+def test_unpack_rejects_odd_vector() -> None:
+    with pytest.raises(ValueError, match="even"):
+        unpack([1.0, 2.0, 3.0])
+
+
+def test_rake_conversion_returns_physical_components() -> None:
+    strike_slip, dip_slip = from_rake([2.0, -3.0], rake_degrees=30.0)
+
+    npt.assert_allclose(strike_slip, np.array([2.0, -3.0]) * np.cos(np.pi / 6))
+    npt.assert_allclose(dip_slip, np.array([2.0, -3.0]) * 0.5)
+    npt.assert_allclose(magnitude(strike_slip, dip_slip), [2.0, 3.0])
+    npt.assert_allclose(rake(strike_slip, dip_slip), [30.0, -150.0])
+
+
+def test_azimuth_conversion_uses_per_patch_strike() -> None:
+    strike_slip, dip_slip = from_azimuth(
+        [1.0, 2.0],
+        azimuth_degrees=90.0,
+        fault_strike_degrees=[0.0, 90.0],
+    )
+
+    npt.assert_allclose(strike_slip, [0.0, 2.0], atol=1e-15)
+    npt.assert_allclose(dip_slip, [1.0, 0.0], atol=1e-15)
+
+
+def test_plate_conversion_round_trip() -> None:
+    strike_slip, dip_slip = from_plate(
+        parallel=[1.0, 2.0],
+        perpendicular=[3.0, 4.0],
+        plate_rake_degrees=[0.0, 90.0],
+    )
+
+    parallel, perpendicular = to_plate(
+        strike_slip,
+        dip_slip,
+        plate_rake_degrees=[0.0, 90.0],
+    )
+
+    npt.assert_allclose(strike_slip, [1.0, -4.0], atol=1e-15)
+    npt.assert_allclose(dip_slip, [3.0, 2.0], atol=1e-15)
+    npt.assert_allclose(parallel, [1.0, 2.0], atol=1e-15)
+    npt.assert_allclose(perpendicular, [3.0, 4.0], atol=1e-15)
+
+
+def test_euler_pole_defines_plate_rake() -> None:
+    fault = Fault.planar(
+        lat=0.0,
+        lon=100.0,
+        depth=10_000.0,
+        strike=0.0,
+        dip=30.0,
+        length=20_000.0,
+        width=10_000.0,
+        n_length=2,
+        n_width=1,
+    )
+
+    plate_rake = plate_rake_from_euler(fault, (90.0, 0.0, 1.0))
+
+    npt.assert_allclose(plate_rake, 90.0, atol=0.2)
diff --git a/tutorials/02_discretization_and_g_matrix.ipynb b/tutorials/02_discretization_and_g_matrix.ipynb
index 986b8c6..a5b50ad 100644
--- a/tutorials/02_discretization_and_g_matrix.ipynb
+++ b/tutorials/02_discretization_and_g_matrix.ipynb
@@ -220,7 +220,7 @@
     "\n",
     "`fault.greens_matrix()` assembles the very same matrix in a single vectorized\n",
     "call — no Python loop over patches. When you build $G$ for a *dataset*,\n",
-    "`geodef.greens.greens(fault, dataset)` does this and then projects each column\n",
+    "`geodef.greens.matrix(fault, dataset)` does this and then projects each column\n",
     "into the data's observation space. For 3-component GNSS that projection is just\n",
     "the `[e, n, u]` interleaving, so the two agree exactly."
    ]
@@ -255,7 +255,7 @@
     "gnss = geodef.GNSS(lon=obs_lon, lat=obs_lat,\n",
     "                   ve=np.zeros(M), vn=np.zeros(M), vu=np.zeros(M),\n",
     "                   se=np.ones(M), sn=np.ones(M), su=np.ones(M))\n",
-    "G_proj = geodef.greens.greens(fault, gnss)\n",
+    "G_proj = geodef.greens.matrix(fault, gnss)\n",
     "print(\"greens(fault, gnss) == greens_matrix:\", np.allclose(G_proj, G))"
    ]
   },
@@ -269,7 +269,7 @@
     "Each entry of $G$ comes from evaluating the Okada solution for one patch–station\n",
     "pair, so assembling $G$ for a large fault and dense data can be slow — and we\n",
     "usually need the *same* $G$ many times (every step of an inversion, every\n",
-    "regularization weight). `geodef.greens.greens()` therefore **hashes its inputs\n",
+    "regularization weight). `geodef.greens.matrix()` therefore **hashes its inputs\n",
     "(fault geometry plus observation points) and caches the result to disk**: the\n",
     "first call computes $G$, and identical later calls just load it. You rarely touch\n",
     "the cache API directly, but it is there if you need it (`geodef.cache.info()`,\n",
@@ -308,9 +308,9 @@
     "cache.set_dir(tempfile.mkdtemp())          # isolated cache dir just for this demo\n",
     "\n",
     "t0 = time.perf_counter()\n",
-    "G_a = geodef.greens.greens(fault, gnss)    # first call: computes G, caches to disk\n",
+    "G_a = geodef.greens.matrix(fault, gnss)    # first call: computes G, caches to disk\n",
     "t1 = time.perf_counter()\n",
-    "G_b = geodef.greens.greens(fault, gnss)    # second call: loaded from disk\n",
+    "G_b = geodef.greens.matrix(fault, gnss)    # second call: loaded from disk\n",
     "t2 = time.perf_counter()\n",
     "\n",
     "print(f\"first call  (compute + cache): {1e3 * (t1 - t0):6.1f} ms\")\n",
@@ -521,7 +521,7 @@
     "## Summary\n",
     "- $G$ is a design matrix: every column is one patch's unit-slip displacement\n",
     "  field.\n",
-    "- `fault.greens_matrix()` (or `geodef.greens.greens()` for a dataset) builds it\n",
+    "- `fault.greens_matrix()` (or `geodef.greens.matrix()` for a dataset) builds it\n",
     "  in one call; we verified it equals the hand-built version.\n",
     "- Columns are blocked `[strike | dip]`; rows are interleaved `[e, n, u]`.\n",
     "- Discretization trades resolution against stability — the motivation for the\n",
diff --git a/tutorials/03_unregularized_inversion.ipynb b/tutorials/03_unregularized_inversion.ipynb
index 8a52e19..02765ba 100644
--- a/tutorials/03_unregularized_inversion.ipynb
+++ b/tutorials/03_unregularized_inversion.ipynb
@@ -14,7 +14,7 @@
     "\n",
     "## Learning objectives\n",
     "- Set up the linear inverse problem $\\mathbf{d} = G\\mathbf{m} + \\boldsymbol{\\varepsilon}$.\n",
-    "- Solve it by (weighted) least squares with `geodef.invert()`.\n",
+    "- Solve it by (weighted) least squares with `geodef.invert.solve()`.\n",
     "- Read the fit with `InversionResult` and `plot.fit`.\n",
     "- See an **unregularized** inversion overfit the noise and oscillate wildly —\n",
     "  the motivation for regularization in notebook 04."
@@ -114,7 +114,7 @@
     "# Synthetic data: forward-model the true slip, then add seeded Gaussian noise.\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
@@ -198,7 +198,7 @@
     }
    ],
    "source": [
-    "raw = geodef.invert(fault, gnss, components=\"dip\")\n",
+    "raw = geodef.invert.solve(fault, gnss, components=\"dip\")\n",
     "print(f\"method:  weighted least squares (unregularized)\")\n",
     "print(f\"reduced chi^2: {raw.reduced_chi2:.3f}\")\n",
     "print(f\"RMS residual:  {raw.rms * 1000:.2f} mm\")"
diff --git a/tutorials/04_regularization.ipynb b/tutorials/04_regularization.ipynb
index 845fadd..da0221c 100644
--- a/tutorials/04_regularization.ipynb
+++ b/tutorials/04_regularization.ipynb
@@ -120,7 +120,7 @@
     "# Synthetic data: forward-model the true slip, then add seeded Gaussian noise.\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
@@ -158,8 +158,8 @@
    },
    "outputs": [],
    "source": [
-    "raw = geodef.invert(fault, gnss, components=\"dip\")                       # lambda = 0\n",
-    "sm = geodef.invert(fault, gnss, components=\"dip\",\n",
+    "raw = geodef.invert.solve(fault, gnss, components=\"dip\")                 # lambda = 0\n",
+    "sm = geodef.invert.solve(fault, gnss, components=\"dip\",\n",
     "                   smoothing=\"laplacian\", smoothing_strength=1.0)\n",
     "\n",
     "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
@@ -226,7 +226,7 @@
     "                 vmin=-vmax, vmax=vmax, colorbar=False, title=\"True\")\n",
     "recovered = []\n",
     "for ax, lam in zip(axes[0, 1:], lambdas):\n",
-    "    r = geodef.invert(fault, gnss, components=\"dip\",\n",
+    "    r = geodef.invert.solve(fault, gnss, components=\"dip\",\n",
     "                      smoothing=\"laplacian\", smoothing_strength=lam)\n",
     "    recovered.append(r.slip_vector)\n",
     "    geodef.plot.slip(fault, r.slip_vector, ax=ax, cmap=\"RdBu_r\",\n",
@@ -284,9 +284,9 @@
     }
    ],
    "source": [
-    "sm = geodef.invert(fault, gnss, components=\"dip\",\n",
+    "sm = geodef.invert.solve(fault, gnss, components=\"dip\",\n",
     "                   smoothing=\"laplacian\", smoothing_strength=1.0)\n",
-    "dm = geodef.invert(fault, gnss, components=\"dip\",\n",
+    "dm = geodef.invert.solve(fault, gnss, components=\"dip\",\n",
     "                   smoothing=\"damping\", smoothing_strength=1.0)\n",
     "\n",
     "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
diff --git a/tutorials/05_choosing_regularization.ipynb b/tutorials/05_choosing_regularization.ipynb
index 9db516a..4694e51 100644
--- a/tutorials/05_choosing_regularization.ipynb
+++ b/tutorials/05_choosing_regularization.ipynb
@@ -103,7 +103,7 @@
     "# Synthetic data: forward-model the true slip, then add seeded Gaussian noise.\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
@@ -277,7 +277,7 @@
     }
    ],
    "source": [
-    "cv = geodef.invert(fault, gnss, components=\"dip\",\n",
+    "cv = geodef.invert.solve(fault, gnss, components=\"dip\",\n",
     "                   smoothing=\"laplacian\", smoothing_strength=\"cv\", cv_folds=5)\n",
     "print(f\"CV-selected:    lambda = {cv.smoothing_strength:.3g}\")"
    ]
@@ -326,7 +326,7 @@
     "                 vmin=-vmax, vmax=vmax, colorbar=False, title=\"True\")\n",
     "recovered = {}\n",
     "for ax, (name, lam) in zip(axes[0, 1:], choices.items()):\n",
-    "    r = geodef.invert(fault, gnss, components=\"dip\",\n",
+    "    r = geodef.invert.solve(fault, gnss, components=\"dip\",\n",
     "                      smoothing=\"laplacian\", smoothing_strength=lam)\n",
     "    recovered[name] = r.slip_vector\n",
     "    geodef.plot.slip(fault, r.slip_vector, ax=ax, cmap=\"RdBu_r\",\n",
diff --git a/tutorials/06_multiple_datasets.ipynb b/tutorials/06_multiple_datasets.ipynb
index 269e7c7..22a8be1 100644
--- a/tutorials/06_multiple_datasets.ipynb
+++ b/tutorials/06_multiple_datasets.ipynb
@@ -15,7 +15,7 @@
     "## Learning objectives\n",
     "- Understand the InSAR **line-of-sight (LOS)** projection and look vectors.\n",
     "- Stack multiple datasets into one linear system.\n",
-    "- Run a **joint** inversion with `geodef.invert(fault, [gnss, insar])`.\n",
+    "- Run a **joint** inversion with `geodef.invert.solve(fault, [gnss, insar])`.\n",
     "- See how relative data weighting shifts the solution."
    ]
   },
@@ -116,7 +116,7 @@
     "# Synthetic data: forward-model the true slip, then add seeded Gaussian noise.\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
@@ -164,7 +164,7 @@
     "look_u = np.full(n_px, 0.92)\n",
     "\n",
     "# Forward-model the true slip into LOS, then add correlated-free noise.\n",
-    "Gi = geodef.greens.greens(\n",
+    "Gi = geodef.greens.matrix(\n",
     "    fault, geodef.InSAR(lon=ilon, lat=ilat, los=np.zeros(n_px), sigma=np.ones(n_px),\n",
     "                        look_e=look_e, look_n=look_n, look_u=look_u)\n",
     ")\n",
@@ -234,9 +234,9 @@
    "outputs": [],
    "source": [
     "kw = dict(components=\"dip\", smoothing=\"laplacian\", smoothing_strength=1.0)\n",
-    "r_gnss = geodef.invert(fault, gnss, **kw)\n",
-    "r_insar = geodef.invert(fault, insar, **kw)\n",
-    "r_joint = geodef.invert(fault, [gnss, insar], **kw)\n",
+    "r_gnss = geodef.invert.solve(fault, gnss, **kw)\n",
+    "r_insar = geodef.invert.solve(fault, insar, **kw)\n",
+    "r_joint = geodef.invert.solve(fault, [gnss, insar], **kw)\n",
     "\n",
     "truth = slip_true[N:]\n",
     "vmax = truth.max()\n",
diff --git a/tutorials/07_correlated_noise.ipynb b/tutorials/07_correlated_noise.ipynb
index dfb7403..37bf2ab 100644
--- a/tutorials/07_correlated_noise.ipynb
+++ b/tutorials/07_correlated_noise.ipynb
@@ -107,7 +107,7 @@
     "# Synthetic data: forward-model the true slip, then add seeded Gaussian noise.\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
@@ -158,7 +158,7 @@
     "n_px = ilon.size\n",
     "look_e, look_n, look_u = np.full(n_px, -0.38), np.full(n_px, 0.09), np.full(n_px, 0.92)\n",
     "\n",
-    "Gi = geodef.greens.greens(\n",
+    "Gi = geodef.greens.matrix(\n",
     "    fault, geodef.InSAR(lon=ilon, lat=ilat, los=np.zeros(n_px), sigma=np.ones(n_px), look_e=look_e, look_n=look_n, look_u=look_u)\n",
     ")\n",
     "\n",
@@ -249,8 +249,8 @@
     "                          covariance=C_full)\n",
     "\n",
     "kw = dict(components=\"dip\", smoothing=\"laplacian\", smoothing_strength=1.0)\n",
-    "r_diag = geodef.invert(fault, insar_diag, **kw)\n",
-    "r_full = geodef.invert(fault, insar_full, **kw)\n",
+    "r_diag = geodef.invert.solve(fault, insar_diag, **kw)\n",
+    "r_full = geodef.invert.solve(fault, insar_full, **kw)\n",
     "\n",
     "fig, axes = plt.subplots(1, 3, figsize=(13, 3.6), constrained_layout=True)\n",
     "vmax = slip_true[N:].max()\n",
diff --git a/tutorials/08_bounds_and_constraints.ipynb b/tutorials/08_bounds_and_constraints.ipynb
index 766dd54..3910bc1 100644
--- a/tutorials/08_bounds_and_constraints.ipynb
+++ b/tutorials/08_bounds_and_constraints.ipynb
@@ -33,7 +33,7 @@
     "- slip magnitude is capped: $\\mathbf m \\le u$;\n",
     "- more general geologic limits: $C\\mathbf m \\le \\mathbf d_{\\text{ineq}}$.\n",
     "\n",
-    "`geodef.invert()` handles these through the `bounds` and `constraints` arguments\n",
+    "`geodef.invert.solve()` handles these through the `bounds` and `constraints` arguments\n",
     "and picks the right solver automatically:\n",
     "\n",
     "- `bounds=(0, None)` → **non-negative least squares (NNLS)**;\n",
@@ -104,7 +104,7 @@
     "# Synthetic data: forward-model the true slip, then add seeded Gaussian noise.\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
@@ -153,8 +153,8 @@
    ],
    "source": [
     "kw = dict(components=\"dip\", smoothing=\"laplacian\", smoothing_strength=0.2)\n",
-    "free = geodef.invert(fault, gnss, **kw)                    # can go negative\n",
-    "nn = geodef.invert(fault, gnss, bounds=(0, None), **kw)    # auto-NNLS\n",
+    "free = geodef.invert.solve(fault, gnss, **kw)                    # can go negative\n",
+    "nn = geodef.invert.solve(fault, gnss, bounds=(0, None), **kw)    # auto-NNLS\n",
     "\n",
     "print(f\"unconstrained min slip: {free.slip_vector.min()*100:+.1f} cm (spurious back-slip)\")\n",
     "print(f\"non-negative  min slip: {nn.slip_vector.min()*100:+.1f} cm\")\n",
@@ -203,7 +203,7 @@
     }
    ],
    "source": [
-    "capped = geodef.invert(fault, gnss, bounds=(0, 1.0), **kw)\n",
+    "capped = geodef.invert.solve(fault, gnss, bounds=(0, 1.0), **kw)\n",
     "print(f\"max slip, uncapped: {nn.slip_vector.max():.2f} m\")\n",
     "print(f\"max slip, capped:   {capped.slip_vector.max():.2f} m (hits the 1.0 m ceiling)\")\n",
     "print(f\"reduced chi^2, uncapped {nn.reduced_chi2:.2f}  vs capped {capped.reduced_chi2:.2f}\")"
@@ -246,7 +246,7 @@
    ],
    "source": [
     "# rake = 90 deg is pure dip-slip in the local frame (matches our true slip).\n",
-    "rk = geodef.invert(fault, gnss, components=\"rake\", rake=90.0,\n",
+    "rk = geodef.invert.solve(fault, gnss, components=\"rake\", rake=90.0,\n",
     "                   smoothing=\"laplacian\", smoothing_strength=1.0)\n",
     "print(f\"components='rake' solves {rk.slip_vector.size} amplitudes \"\n",
     "      f\"(vs {2*N} for 'both')\")\n",
diff --git a/tutorials/09_uncertainty_and_resolution.ipynb b/tutorials/09_uncertainty_and_resolution.ipynb
index 357cab8..0ec1742 100644
--- a/tutorials/09_uncertainty_and_resolution.ipynb
+++ b/tutorials/09_uncertainty_and_resolution.ipynb
@@ -107,7 +107,7 @@
     "# Synthetic data: forward-model the true slip, then add seeded Gaussian noise.\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
@@ -134,7 +134,7 @@
    },
    "outputs": [],
    "source": [
-    "result = geodef.invert(fault, gnss, components=\"dip\",\n",
+    "result = geodef.invert.solve(fault, gnss, components=\"dip\",\n",
     "                       smoothing=\"laplacian\", smoothing_strength=1.0)"
    ]
   },
@@ -243,7 +243,7 @@
     "d_check = G_full @ m_check + rng.normal(0.0, sigma, G_full.shape[0])\n",
     "gnss_check = geodef.GNSS(lon=glon, lat=glat, ve=d_check[0::3], vn=d_check[1::3], vu=d_check[2::3],\n",
     "                         se=np.full(n_sta, sigma), sn=np.full(n_sta, sigma), su=np.full(n_sta, sigma))\n",
-    "rec = geodef.invert(fault, gnss_check, components=\"dip\",\n",
+    "rec = geodef.invert.solve(fault, gnss_check, components=\"dip\",\n",
     "                    smoothing=\"laplacian\", smoothing_strength=0.3)\n",
     "\n",
     "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))\n",
diff --git a/tutorials/10_nonlinear_geometry.ipynb b/tutorials/10_nonlinear_geometry.ipynb
index b42fc67..45df6b6 100644
--- a/tutorials/10_nonlinear_geometry.ipynb
+++ b/tutorials/10_nonlinear_geometry.ipynb
@@ -111,7 +111,7 @@
     "# Synthetic data: forward-model the true slip, then add seeded Gaussian noise.\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
@@ -155,7 +155,7 @@
     "    )\n",
     "\n",
     "def misfit(dip):\n",
-    "    r = geodef.invert(make_fault(dip), gnss, components=\"dip\",\n",
+    "    r = geodef.invert.solve(make_fault(dip), gnss, components=\"dip\",\n",
     "                      smoothing=\"laplacian\", smoothing_strength=1.0)\n",
     "    return r.reduced_chi2"
    ]
@@ -274,7 +274,7 @@
    "outputs": [],
    "source": [
     "best_fault = make_fault(opt.x)\n",
-    "best = geodef.invert(best_fault, gnss, components=\"dip\",\n",
+    "best = geodef.invert.solve(best_fault, gnss, components=\"dip\",\n",
     "                     smoothing=\"laplacian\", smoothing_strength=1.0)\n",
     "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))\n",
     "vmax = slip_true[N:].max()\n",
diff --git a/tutorials/11_gradient_geometry.ipynb b/tutorials/11_gradient_geometry.ipynb
index ba49d33..3675140 100644
--- a/tutorials/11_gradient_geometry.ipynb
+++ b/tutorials/11_gradient_geometry.ipynb
@@ -64,7 +64,7 @@
     "\n",
     "_zero = np.zeros(n_sta)\n",
     "_one = np.ones(n_sta)\n",
-    "G_full = geodef.greens.greens(\n",
+    "G_full = geodef.greens.matrix(\n",
     "    fault, geodef.GNSS(lon=glon, lat=glat, ve=_zero, vn=_zero, vu=_zero, se=_one, sn=_one, su=_one)\n",
     ")\n",
     "sigma = 0.01  # 1 cm station noise\n",
diff --git a/tutorials/OUTLINE.md b/tutorials/OUTLINE.md
index f7fabf1..308c5ae 100644
--- a/tutorials/OUTLINE.md
+++ b/tutorials/OUTLINE.md
@@ -118,7 +118,7 @@ So notebooks stay consistent and composable:
   `np.allclose(...)`. This demystifies the API without hiding it. Reserve it for
   the rare spots where the manual version is itself instructive — e.g. `d = G @ m`
   vs. `fault.displacement()` in Tutorial 01, or building `G` column by column vs.
-  `greens.greens()` in Tutorial 02. Skip it everywhere else: the default remains
+  `greens.matrix()` in Tutorial 02. Skip it everywhere else: the default remains
   one clear `geodef.*` call plus a labeled plot.
 
 ---
@@ -214,7 +214,7 @@ structure as the discrete forward operator.
 - How dataset projection turns the raw 3-component response into observed rows
   (interleaved E/N/U for GNSS; foreshadow LOS for InSAR in Tut 06).
 - **Double-demo:** build `G` column by column (unit slip per patch) and confirm
-  it equals `fault.greens_matrix(...)` / `geodef.greens.greens(...)`.
+  it equals `fault.greens_matrix(...)` / `geodef.greens.matrix(...)`.
 - **Sidebar (absorbs old caching notebook):** assembling `G` is expensive and
   often repeated, so GeoDef hashes its inputs and caches `G` to disk; one short
   timing demo (first call computes, second loads).
@@ -222,7 +222,7 @@ structure as the discrete forward operator.
   teaser for Tutorials 03–04.
 
 **Key calls.** `Fault.planar(...)` with multiple patches, `fault.greens_matrix`
-and `geodef.greens.greens(fault, dataset)` (shown to agree with a hand-built
+and `geodef.greens.matrix(fault, dataset)` (shown to agree with a hand-built
 `G`), `np.linalg.lstsq` for the warm-up, `fault.patch_index`, the
 `geodef.cache.info()` / `set_dir()` sidebar, `plot.slip`, `plot.vectors`.
 
@@ -255,7 +255,7 @@ inversion overfit noise.
   watch the recovered slip oscillate wildly while fitting the data "too well."
   Motivates Tutorial 04.
 
-**Key calls.** `geodef.invert(fault, dataset)` (default WLS),
+**Key calls.** `geodef.invert.solve(fault, dataset)` (default WLS),
 `InversionResult` (`.slip_vector`, `.predicted`, misfit fields), `plot.fit`,
 `plot.slip`.
 
@@ -287,7 +287,7 @@ common regularization operators.
   "equations" with weight `√λ`.
 - Qualitative effect of `λ`: under- vs. over-smoothing (sets up Tutorial 05).
 
-**Key calls.** `geodef.invert(..., smoothing='laplacian'|'damping'|'stresskernel',
+**Key calls.** `geodef.invert.solve(..., smoothing='laplacian'|'damping'|'stresskernel',
 smoothing_strength=λ, smoothing_target=m_ref)`; `greens` Laplacian builder
 referenced conceptually.
 
@@ -316,7 +316,7 @@ guess a good value (then check it in Tutorial 05).
 - When the methods agree/disagree and how to choose between them.
 
 **Key calls.** `geodef.lcurve(...)`, `geodef.abic_curve(...)`,
-`geodef.compute_abic(...)`, and `geodef.invert(..., smoothing_strength='abic'`
+`geodef.compute_abic(...)`, and `geodef.invert.solve(..., smoothing_strength='abic'`
 `|'cv', cv_folds=...)`; their built-in curve plots.
 
 **Plots.** L-curve with marked corner; ABIC vs. `λ`; CV error vs. `λ`; the
@@ -344,7 +344,7 @@ relative weighting.
   connection to the regularization hyperparameter (foreshadow multi-λ).
 
 **Key calls.** `geodef.GNSS(...)`, `geodef.InSAR(...)` with look vectors,
-`geodef.invert(fault, [gnss, insar], ...)`, `plot.insar`, `plot.vectors`.
+`geodef.invert.solve(fault, [gnss, insar], ...)`, `plot.insar`, `plot.vectors`.
 
 **Plots.** GNSS vectors and InSAR LOS for the same scenario; joint-inversion
 slip vs. single-dataset slip; per-dataset fit panels.
@@ -375,7 +375,7 @@ InSAR.
 - Practical down-sampling of dense InSAR as a related concern (brief).
 
 **Key calls.** `geodef.InSAR(...)` with a full covariance / covariance-function
-specification; `geodef.invert(...)` with the resulting `C_d`.
+specification; `geodef.invert.solve(...)` with the resulting `C_d`.
 
 **Plots.** A covariance matrix / covariance function; inversion with diagonal
 vs. full `C_d` and the difference in recovered slip and uncertainty.
@@ -402,7 +402,7 @@ slip direction.
   (`components='azimuth'`), which also encodes a sign/sense prior cleanly.
 - Trade-offs: constraints vs. smoothing; bias vs. admissibility.
 
-**Key calls.** `geodef.invert(..., bounds=(0, None))` (auto-NNLS),
+**Key calls.** `geodef.invert.solve(..., bounds=(0, None))` (auto-NNLS),
 `bounds=(lb, ub)` (bounded LS), `method='constrained', constraints=(C, d)`,
 `components='rake', rake=...`, `components='azimuth', slip_azimuth=...`.
 
@@ -465,7 +465,7 @@ on top of the linear slip inversion.
   (`emcee`) for posterior uncertainty on geometry — pointer to a future
   `examples/` study.
 
-**Key calls.** A small Python objective wrapping `geodef.invert(...)` inside
+**Key calls.** A small Python objective wrapping `geodef.invert.solve(...)` inside
 `scipy.optimize.minimize` / a grid loop; reuse of earlier inversion calls.
 
 **Plots.** Misfit vs. a scanned geometry parameter (e.g. dip); recovered vs.
diff --git a/tutorials/reference_plots.ipynb b/tutorials/reference_plots.ipynb
index 7ad3cac..aa60490 100644
--- a/tutorials/reference_plots.ipynb
+++ b/tutorials/reference_plots.ipynb
@@ -104,7 +104,7 @@
     ")\n",
     "\n",
     "# --- Inversion ---\n",
-    "result = geodef.invert(\n",
+    "result = geodef.invert.solve(\n",
     "    fault, [gnss, insar],\n",
     "    smoothing=\"laplacian\", smoothing_strength=1e3,\n",
     ")\n",