An n-simplex is the convex hull of n+1 points that do not lie in some sub-space of less than n dimensions. I use the term zone for a n-simplex in n dimensional space. A simplicial complex is a partition of space into zones. Poincaré introduced the simplicial complex in order to reason about topology. Tulio Regge proposed partitioning a Riemann space into simplex shaped zones and as an approximation, calling each zone flat, i.e. Euclidean. If the length of each edge is specified then the shape of each zone is determined and thus also the dihedral angles. The dihedral angle is the angle between two adjacent faces of a zone.

A 3D simplex is a tetrahedron. The dihedral angle of a triangle is the ordinary angle at a vertex.

Regge pointed out that if you sweep up all of the curvature of n-space into the smallest possible piles, the piles are n−2 dimensional. He called these bones of curvature. For his partitioning into zones, the bones are (n−2)-simplexes. The curvature of the surface of a 3D polyhedron is at the vertices.

The curvature in our Riemannian space can be computed by summing the dihedral angles about a n−2 dimensional bone and subtracting 2 π. Thus the curvature at the corner of a cube is 2 π − 3*(right angle of a square) = π/2.

One necessary elementary geometry exercise is to compute a geodesic. Within a zone it is an ordinary straight line since the zone is flat. The trick is to carry the line into the neighboring zone.

If we consider two zones sharing a (n−1)-simplex as a boundary then each will have a unique unshared vertex that we will take to be the origin of a coordinate system.
We take the n edges emanating from the origin to be a bases for a coordinate system for that zone.
First we observe that if we know the square of the lengths of the edges then it is trivial to compute the metric tensor for this coordinate system.
g_{ij} = edge_{i} • edge_{j} .
The dot product of sides a and b of a triangle is (a^{2} + b^{2} − c^{2})/2.

We can embed two neighboring zones in one flat n dimensional space and compute the transformation between these two coordinate systems. This will allow us to transform vectors and tensors at the interface from one coordinate system to the other. It seems strategic to compute the coordinates of second origin with respect to the first origin.

Let x^{i} be the coordinates of the second origin with respect to the first coordinate system.
Let a_{i} be the known lengths of the basis vectors of the second coordinate system.
Referring to points by means of their coordinates in the first system, basis vector k of the second system is the line segment from x^{i} to δ^{ik}.
Its components are (δ^{i}_{k} − x^{i}).
Its squared length is g_{ij}(δ^{i}_{k} − x^{i})(δ^{j}_{k} − x^{j}).
g_{ij} is the metric tensor of the first system in this very funny looking tensor expression.
The summation is over i and j.
δ^{ij} is the Kronecker delta and provides the coordinates of the vertices of the shared simplex that is the boundary between our two zones.
We must thus solve for each x^{i} in the n equations:

g_{ij}(δ^{ik} − x^{i})(δ^{jk} − x^{j}) = a^{k}

g_{ij}δ^{ik}δ^{jk}
− g_{ij}δ^{ik}x^{j}
− g_{ij}x^{i}δ^{jk}
+ g_{ij}x^{i}x^{j} = a^{k}

(g_{kj}+g_{jk})x^{j} − g_{ij}x^{i}x^{j}
= g_{kk} − a^{k}

2g_{kj}x^{j} − g_{ij}x^{i}x^{j}
= g_{kk} − a^{k}

Recall that the g’s are all known.
Subtracting these equations from each other for distinct values of k provides n−1 independent linear equations which we can use to eliminate all but one x.
What remains is a quadratic equation in that x.
The two solutions lead to locating the secondary origin at two points equidistant from the interface between the zones.
The sum of the x^{i} will exceed 1 one for the correct choice of the root and be less by the same amount for wrong choice.

With the coordinates of the second origin in terms of the first, it is easy to compute the linear transformation between the coordinate systems.
If we circulate around a curvature bone, multiplying by this coordinate transformation each time we go from one zone to the next, we arrive back where we started with a composed transformation that will be a simple rotation about the bone.
If that rotation is w^{ij}, then the summation of w^{ij}w^{kl}, over the five zone faces (each a bone) will provide an approximation to R^{ijkl}, the general Riemann curvature tensor.
(I think that it is necessary to weight the sum by the 3D volume of the face.)
I have glossed over the fact that there are several coordinate systems for each zone, one with each of the five vertices serving as origin, but the transformation between these is immediate.

Note that we can do ordinary special relativistic physics at this point, for we can compare and transport vectors and tensors across zone boundaries. I have not yet described how to impose the equations of General Relativity. That will necessarily require governing the edge lengths which are viewed as inputs to the calculations suggested above. I am inclined to write code to do the above calculations for it will serve to debug code that performs the initial value problems of GR.

Here are more concrete notes on code to do this. Here is some Mathematica code that could be used to debug the calculations suggested above. There is a great deal of literature on the web on this subject that can be found by looking for “Regge calculus” thru Alta Vista or Google.