ymf mesh tags and provenance, coefficients that name unknowns, domains beyond a box (cekees/ymf#6-#8)

docs/README.md

@@ -21,6 +21,7 @@ and lands in `design/` with the code that implements it.
 | Note | Status |
 |---|---|
 | [design/time-outside-space.md](design/time-outside-space.md) | **Implemented** 2026-10-26: archive format 2, one list of instants with the spaces inside each. |
+| [design/domains-beyond-box.md](design/domains-beyond-box.md) | **Implemented** 2026-11-04 (C27 Phase 5): boundary regions by tag, periodic on a region whose faces match, `file` and `discrete` domains; decisions D1-D4. |
 
 ## History
 

docs/archive-format.md

@@ -398,6 +398,89 @@ per-rank grid metadata as these strings and stores them in fixed-width HDF5
 datasets. (Parallel HDF5 rejects variable-length strings, so the width is
 agreed across ranks with an `allreduce` of the maximum length.)
 
+## Meshes: facets, provenance and the cell type's order
+
+A mesh is a grid: a document with one instant whose one space is the mesh
+(proteus names it `Mesh Spatial_Domain`, as in its run archives). Two
+optional keys, on a uniform space or on any grid of a spatial collection,
+carry what only a mesh needs:
+
+```yaml
+topology: {type: Triangle, number_of_elements: 2, data_item: ...}
+geometry: {type: XYZ, data_item: ...}
+attributes: [...]
+facets:
+  topology: {type: Polyline, number_of_elements: 4, nodes_per_element: 2,
+             data_item: ...}          # node numbers of each facet
+  tags: {data_item: ...}              # one integer per facet
+  view: global                        # or subdomain
+provenance:
+  tool: a mesh converter              # free text: the archive names no tools
+  version: '0.1'
+  files: [{name: square.msh, sha256: 3f5a...}]
+  extra: {...}                        # anything else the writer keeps
+```
+
+**Facets** are the tagged faces of the mesh (edges in 2D, points in 1D),
+listed by their node numbers, with an integer tag each: **0 untagged, >0 a
+condition**. The names of the conditions belong to the specification; the
+mesh holds integers. A tag per facet, not per node, because a node at a
+corner belongs to two faces that may carry different conditions. The facet
+type follows the cell type and has the same order: `Polyline` facets for
+`Triangle` and `Quadrilateral` cells, `Edge_3` for `Tri_6`, `Triangle` for
+`Tetrahedron`, `Tri_6` for `Tet_10`, `Quadrilateral` for `Hexahedron`
+(wedges and pyramids may use `Mixed`).
+**-1 marks a face on the interface between two ranks.** It is the
+partitioner's product, so a mesh as a whole (`view: global`, the default)
+never stores it; one rank's piece (`view: subdomain`) may. `write_ymf`
+refuses a -1 in a global view, and any tag below -1, when the tags are
+inline; `ymf.mesh.check_facets(grid, base_dir)` reads tags held in an HDF5
+file and also checks that every facet node is a node of the grid and every
+facet a face of some cell.
+
+**Provenance** says where the mesh came from: the tool, its version, and
+each source file's name and SHA-256 (`ymf.mesh.file_provenance(tool, paths)`
+computes them), so a reader can tell whether the mesh still matches the
+files beside it. It takes no other keys; `extra` holds the rest.
+
+**Geometry order is the cell type's.** `Tri_6` is a quadratic triangle and
+`Tet_10` a quadratic tetrahedron; there is no separate order field to
+disagree with the type (`geometry_order(type)` reads it off). A quadratic
+cell lists its vertices first, in the linear type's order, then one node on
+each edge, in the order of `QUADRATIC_EDGES`:
+
+| type | edge nodes after the vertices sit on |
+|---|---|
+| `Edge_3` | (0,1) |
+| `Tri_6` | (0,1) (1,2) (2,0) |
+| `Quad_8` | (0,1) (1,2) (2,3) (3,0) |
+| `Tet_10` | (0,1) (1,2) (2,0) (0,3) (1,3) (2,3) |
+| `Pyramid_13` | (0,1) (1,2) (2,3) (3,0) (0,4) (1,4) (2,4) (3,4) |
+| `Wedge_15` | (0,1) (1,2) (2,0) (3,4) (4,5) (5,3) (0,3) (1,4) (2,5) |
+| `Hex_20` | (0,1) (1,2) (2,3) (3,0) (4,5) (5,6) (6,7) (7,4) (0,4) (1,5) (2,6) (3,7) |
+
+This is VTK's order, which is what a viewer draws: ParaView's XDMF readers
+pass the connectivity through unpermuted. Sources do not all agree with it.
+A common mesh generator numbers the last two edge nodes of a 10-node
+tetrahedron (2,3) then (1,3); proteus's own quadratic tetrahedron numbers
+them (0,1) (1,2) (2,3) (0,2) (1,3) (0,3). A writer states its source's
+convention as data and permutes:
+
+```python
+from ymf.mesh import reorder_cells, edge_midpoint_offsets
+source = ((0, 1), (1, 2), (2, 0), (0, 3), (2, 3), (1, 3))
+cells = reorder_cells(cells, "Tet_10", source)          # columns 8 and 9 swap
+assert edge_midpoint_offsets("Tet_10", nodes, cells).max() < 1e-12   # straight-sided sample
+```
+
+`edge_midpoint_offsets` is the check: on a straight-sided mesh every edge
+node sits at its edge's midpoint once the order is right, and about half an
+edge away when two are swapped.
+
+`ymf2xmf` writes the grid as before; the facet list and the provenance are
+not written to the `.xmf` (XDMF's `Set` is not supported, see
+`docs/xdmf-model.yaml`).
+
 ## What the writer checks, and what it doesn't
 
 YMF splits checking across a trust boundary:

docs/changes.md

@@ -7,6 +7,165 @@ rediscover why a thing is the way it is. Open gaps are in
 
 ## 2026-11
 
+- **2026-11-04: `discrete` domains are posed (D4).** A mesh as ymf stores
+  one (`topology`, `geometry`, and the tagged `facets` and `provenance` of
+  cekees/ymf#6) is a domain: `adr_problem` poses it whole, with bounds from
+  its inline nodes (or `bounds`, or the solver's, for nodes in a file), its
+  boundary regions by facet tag and `mesh: {levels}` only; periodic is
+  refused. The design note, decided 2026-11-04, lands as
+  `docs/design/domains-beyond-box.md`.
+- **2026-11-04: `file` domains are posed (D2, D3).** `adr_problem` carries
+  a file domain whole (`domain.region`) for the solver to read. Optional
+  `bounds: {lower, upper}` set the box; without them `lower` and `upper`
+  are None and the solver takes them from the mesh (posing follows meshing
+  for this kind alone), and the dimension must be given. A mesh file
+  (`2dm`, `3dm`, `hex`, `tetgen`, `ymf`) fixes the mesh: its discretization
+  gives `mesh: {levels}` only, keyed `levels=L`, and a mesh size is
+  reported as lacking; `gmsh_geo` is meshed by `h`. `ymf.discretizations`,
+  `ymf.closure.realizations` and `ymf.state` take the domain for this
+  (optional; without it a mesh needs a size, as before). Periodic
+  conditions on a file domain are refused.
+- **2026-11-04: periodic conditions on a `piecewise_linear` region (D3).**
+  Refused before ("periodic faces are a box's"). Allowed now when, on each
+  periodic axis, the region's boundary segments (2-D) or facets (3-D) on the
+  lower face translate exactly onto those on the upper face (as point sets,
+  to 1e-8 of the box); otherwise refused with the count of pieces without a
+  partner. `ymf.symbolic.problem.periodic_faces_match`, used by both
+  readers. The solver's part (no boundary points added by the mesher, a
+  partner check after meshing) is Proteus's companion PR.
+- **2026-11-04: boundary regions by tag (C27 Phase 5, decision D1).**
+  `boundary_regions` entries take `tags: [ints or names]` instead of
+  `where`, exactly one of the two; names resolve through the domain's
+  `boundary_tags`. `adr_problem` and `split_problem` give a tag region as
+  `flags` and no `where`. Checked on load: one of the two, every name
+  resolves, no tags on a box, a periodic region uses `where`.
+- **2026-11-04: the shallow-water family.** `tests/fixtures/shallow_water`:
+  the model (depth and discharge; porosity, Manning friction and
+  Darcy-Forchheimer drag as coefficients of the state; `wall` and
+  `outflow` condition sets), Proteus's numerics (SW2DCV as one sub-problem,
+  its algebraic scheme in the free-form `discretization` block), a lake at
+  rest over a bump, and the Torino flume (dam break down a 3 % slope into a
+  rod canopy, from `torino_flume/dswe/torino.py`). The second family
+  needed no new schema: what two-phase flow added covers it.
+- **2026-11-04: boundary conditions in `split_problem`.** Each sub-problem
+  carries the conditions it imposes, in `adr_problem`'s plain-data shape
+  (`where` code, per-component values, series, robin parts, periodic axes),
+  for every kind: those on its own unknowns, and those on a coefficient its
+  introduced unknown carries. In the dam break: flow gets the free-slip
+  walls, `p = 0` and an open `v` on top; VOF gets air flowing in through
+  `χh`; the level-set solves get none, as Proteus's `setAtmosphere` sets
+  none for them.
+- **2026-11-03: gauges.** `model.gauges` in the shape the force-history
+  plan (C22, step 4) set out: `point`, `line` and `line_integral` gauges of
+  named fields, and `boundary_integral` gauges of wetted area, forces and
+  moments per boundary region, with `sample_rate` and `file`. What each kind
+  needs is checked on load. Also: a reference error now says "an analytical
+  solution", not "a analytical solution".
+- **2026-11-03: condition sets.** `strong_form.condition_sets` names bundles
+  of conditions without a region (a model's `free_slip`, `atmosphere`), and
+  `strong_form.boundary_sets` applies them to regions. They are a separate
+  key rather than a `set:` entry among `boundary_conditions` because
+  strictyaml cannot accept either of two mapping shapes in one list, and
+  making a condition's keys optional would lose the schema error that names
+  the file holding an incomplete one. `ymf.boundary_conditions` expands
+  them, and `adr_problem`, `ymf.units` and the load-time checks all read
+  through it; a condition written out replaces a set's member for the same
+  region and variable.
+- **2026-11-03: `slip`, `open`, `inflow`, `characteristic`; conditions on
+  a coefficient.** The boundary kinds two-phase flow and shallow water need:
+  a free-slip wall, an open boundary (C12's declared open boundary), a value
+  imposed where the flow enters, a stage imposed through the incoming
+  characteristics. A condition may name a coefficient of the state (`χ`).
+  Conditions are now checked on load: the region exists, the variable is an
+  unknown or a coefficient, `slip` is on a vector, and a value is given
+  exactly when the kind takes one. `ymf.units` checks `slip`, `inflow` and
+  `characteristic` like `dirichlet`. `adr_problem` refuses the four kinds
+  and conditions on coefficients, naming the solver they are for.
+- **2026-11-03: each sub-problem classified on its own
+  (`ymf.symbolic.split.split_problem`).** The plain data of a split
+  discretization, one ADR form per sub-problem. Its own unknowns are
+  unknowns; what it reads from the others is a known field (a function of
+  the coordinates and time, written by name, its gradient
+  `grad_<field>_<axis>`), and must be declared in `lags`; a coefficient
+  lagged via the sub-problem that carries it is replaced by the carrier.
+  `binds` are substituted (a symbol with only a meaning, `h`, stays a
+  parameter the solver supplies), and `regularize` smooths `heaviside` and
+  `sign` exactly as Proteus does (`gf.H`, `gf.D`, `2 H_ε − 1`). Schema:
+  a sub-problem's `for` and `regularize`, a lag's `of`. The two-phase dam
+  break classifies in all five sub-problems: flow as variable-density
+  Navier-Stokes with ρ and μ from VOF's carrier, VOF and the level set as
+  transport in the lagged velocity (`v·∇φ` is a Hamiltonian in `φ` alone
+  there), redistancing as an eikonal Hamiltonian, mass correction as a
+  regularized reaction-diffusion.
+- **2026-11-03: coefficients of the state keep their names.**
+  `adr_problem` lists each coefficient that follows the state under
+  `coefficients`: its numpy code, `depends_on`, and its derivative in each
+  unknown it depends on (`sharp: true` instead, for `heaviside` or `sign`
+  of an unknown). The ADR form keeps them substituted, so a solver of one
+  system is unaffected; the names are for a generator or a dispatcher. A
+  problem without such coefficients has no `coefficients` key, so its
+  record is unchanged.
+- **2026-11-03: a discretization may be split into sub-problems.**
+  `subproblems:` (solved in order each step: what each `solves`,
+  `enforces`, `introduces`, `lags`, `writes_back`, with its own space,
+  time integration, stabilization and a `solver` that refines the
+  entry's), `binds:` (symbols such as `h` and `ε = 1.5 h` that the
+  discretization supplies), `regularize:` (sharp functions made smooth) and
+  `step_control:`. The structure and every name are checked on load
+  (`ymf.references`); `ymf.discretizations.subproblem_solver` gives a
+  sub-problem's effective solver. Classification per sub-problem is not
+  done yet. Splitting is not a new `method`: a split discretization is
+  still a mesh one, and its studies and output keys work as before.
+- **2026-11-03: `heaviside`, `sign`, `min`, `max`, `norm`, `sym`.** The
+  two-phase model writes `χ = heaviside(φ)`, its initial signed distance
+  with `min`, `max` and `sqrt`, its viscous stress with `sym(∇v)`;
+  shallow-water friction uses `norm(q)`. They parse, print back
+  (`to_notation`), and become numpy code (`min`/`max` as nested
+  `numpy.minimum`/`numpy.maximum`, not `functools.reduce`). The derivative
+  of `heaviside` or `sign` is a Dirac delta, which has no value at a point:
+  where the ADR form needs one (a sharp function of an unknown), it is
+  refused with a message saying the discretization must regularize it.
+- **2026-11-03: named equations.** `strong_form.equations` takes
+  `{name, for, equation}` beside bare strings. Bare strings are still
+  ordered like the unknowns; with `for` on every equation the order written
+  stops mattering. Names merge across `extends` like every other named
+  list, and are what a discretization's sub-problems will cite. Checked on
+  load: a name used twice, `for` naming something that is not an unknown,
+  two equations for one unknown, `for` on some but not all, and (outside an
+  abstract file) an unknown no equation is for. `ymf.equations_in_order`
+  is the one reader, shared by `ymf.symbolic` and `ymf.units`.
+- **2026-11-03: a coefficient's formula may name the unknowns.** A
+  coefficient could be a formula in the coordinates, time and the other
+  coefficients; one that named an unknown was "could not be resolved". Now
+  it may: a state-dependent density, conductivity or friction law is a
+  coefficient whose formula names the unknowns, substituted into the
+  equations so the ADR form carries the dependence and its derivatives.
+  Decided by cekees for the model-family work (two-phase level-set flow:
+  `χ = heaviside(φ)`, `ρ(χ)`, `μ(χ)`; shallow water: Manning friction,
+  canopy drag), in place of a separate `relations:` block: one kind of
+  thing, a coefficient, whether it is constant, varies in space or follows
+  the state. What came with it: boundary data, initial conditions and exact
+  solutions that lean on such a coefficient are refused (they are evaluated
+  at points in x, y, z, t); `ymf.units` warns when a characteristic scale or
+  a dimensionless number is built on one; a coefficient may not take an
+  unknown's name; a formula may not take a derivative of an unknown (yet).
+  `ymf.symbolic.state_coefficients` lists which coefficients follow the
+  state, and on what.
+- **2026-11-03: a mesh carries its tagged facets, its provenance, and its
+  order in its cell type.** Proteus's mesh campaign (C7.1, step A) needs a
+  mesh to say which faces carry which boundary condition, and where it came
+  from, so converters from other codes' inputs can write ymf. A grid may
+  now carry `facets` (node tuples and one integer tag per facet: 0
+  untagged, >0 a condition, -1 a rank interface, allowed only in a
+  `subdomain` view) and `provenance` (tool, version, files with SHA-256;
+  no list of tools anywhere). Facets, not node tags, because a corner node
+  belongs to two faces. Geometry order stays the topology type's
+  (`geometry_order`), and the node order of each quadratic type is now
+  written down (`QUADRATIC_EDGES`, VTK's, what viewers draw); the
+  10-node tetrahedra of one widely used generator and of proteus number
+  their edge nodes differently, and `ymf.mesh.edge_node_permutation`
+  maps any stated convention onto ymf's. Older documents are unchanged.
+
 - **2026-11-01: one name per thing: `model` and `approximation`.** The spec's
   top-level key has been `model` since a66ea1c, and the archive of a spec
   stores the realized spec under the key the spec uses, so it has written