This is to explain that code. The code has the limited goal of supporting ray tracing. This requires the numerical support for an affine connection which is knowing the coordinate transformation across facets.

The Zone module

Routine Zone.newz is called once per zone, during morphogenesis. Arguments to Zone.newz: The mutable neigh has not generally acquired all values upon entry to newz but morphogenesis is the convenient time to deliver access to the array. Presumably the array will be completed before ray is called, and constant thereafter. There are always safe values there to call. An access value is a record of routines to call a neighboring zone depending on the purpose of the message. In current version it is the record:
ray
Introduce geodesic
curl
Curvature query concerning shared bone.

“zbcc” is a documentation term referring to barycentric coordinates in a zone. “fbcc” is a documentation term referring to barycentric coordinates in an n−1 facet, with basis including the natural normal.

Zone.newz computes the covariant, gd, and contravariant, gu, metric tensors and keeps those within the zone. gu is the n+1 by n+1 matrix; each row and column sum to 0. The contravariant metric tensor for choice of any other vertex as origin may be had by eliminating some row and column of gu. ufn is the array of n+1 normals to the facets in zbcc.

Messages between neighboring zones use the numbering plan of their shared facet; the zones do not have access to the global numbering plan. The numbering plan of a zone’s facet is like that of the zone except for eliding the number by which the zone knows the facet. Zone.newz is of type
float array array -> fray nAcc array -> (unit -> 'a) -> int -> fray nAcc.

Values established in the zone during morphogenesis:

Unitized facet normals don’t generally sum to 0!!!!

The real contravariant metric tensor, gij = gu.(i).(j). Omitting row j and column j from gu yields the contravariant metric tensor as if vertex j were chosen as origin.

The arguments to rayIn are:

fn : int
The number of the facet thru which the ray enters.
lam : float
The current value of λ, the geodesic parameter. λ increases from zone to zone. Perhaps the ray will stop when it reaches zero.
wh : float array (facet bcc, (n−1))
Where on the facet the ray enters.
vel : float array (facet bcc, (n), contravariant)
Velocity of ray.
rayIn returns the egress pair in global coordinates. To follow the plan described here to express vectors in facet coordinates we must formalize that concept.

I find the two trivial routines ex and com confusing. I have had more bugs in such logic than all of the linear algebra.

let ex f n = if f>n then n else n+1
When we have the number, n, of a vertex in the numbering plan of facet f, then ex f n yields the number of the same vertex in the zone’s numbering plan.
let cmp f z = if f=z then raise Hit4; if f<z then f else f-1
When we have a local vertex (or facet) number z in the zone numbering plan, and need the number for the same vertex in the numbering plan of facet f, then write cmp f z.

Curvature of Complex

Curvature is measured about a bone which is one of the n−2 simplexes. Such bones are known to the zone by a pair of local vertex numbers. In a zone, any two distinct local vertex indices denote a bone and each bone is denoted by both such a pair, and that pair reversed. Each zone shares n(n+1)/2 bones with other zones. The two facets named by these numbers share the bone as a sub-simplex. The idea is to send a message across one of these facets to the neighbor which indeed shares the bone. The message specifies which of the n bones that are shared between these two neighbors the query concerns. This specification is in the numbering plan of the shared facet. Also in the message are n, or n+1 (design decision) vectors to be transported around the bone by ‘parallel transport’. They are transmitted in zbcc across a facet. When the get back to the originator the message is recognized and the scalar deficit is computed and returned. Actually I think it is a tensor deficit; we shall see. The various zones are free to copy this returned value and pass the result back to the original caller.

In curl code k is the zone’s number for the vertex opposite to which the message will be relayed. That number together with fn which identifies the neighbor who sent the message, constitute the two numbers that identify the shared bone. Expression (if fn < k then fn else fn-1) computes the vertex number, in the facet’s NP, of the facet vertex opposite the bone.

Some Ocaml Types

   type fray = float array * float array
   type 'a rayp = float -> fray -> 'a
   type curlp = int * float array array -> float array array
   type 'a nAcc = {ray : 'a Zone.rayp; curl : Zone.curlp}
fray is a pair of vectors in fbcc that describe the position and direction of a ray or geodesic. rayp is the type of the function called by zone code to introduce a ray to a neighboring zone. The type parameter 'a is the type of the value returned which is of no concern to the zone logic. curlp is the type of function called by a zone to send a message to a neighboring zone on a trip about a bone.

Discarded Ideas

The above was once written around the idea that during morphogenesis each facet would develop a linear transformation between the bcc of the two adjacent zones. I think now this is a bad idea. It is bad because such a transformation is the product of two transformations between zone and facet bccs and each of these are generally sparse. It may be both conceptually simpler and faster to perform the two transformations within the respective zones. The facets would be limited to flipping the sign of the normal and permuting the vertex numbers. This may save time, space and hair.