This is a reprise of code written early in 1956 for an IBM 650 to explore the equations for Lagrangian hydrodynamics with triangular zones in what is now called ‘finite element analysis’. It was the first triangular zoning at Livermore and quite possibly the first anywhere. It was taken then to be a variation of difference equations. The data were scaled fixed point then, but floating point today. Array references were via self modifying code; there were no index registers. The physics is isothermal hydrodynamics.

The mesh size was 20 by 20 which with the program fit in the 2000 words of drum memory.

The type pt holds the state of a pair of triangular zones. x, y are the coordinates of the lower left of the first zone. xd, yd are the velocity components of that point.

Then, as now, we wrote x and y components separately doubling much of the code.

To mostly debug the code we adopt a symmetric initial boundary conditions and verify that the symmetry persists.

dx, dy
Now as then we trade a subtract for a load and store from an array of vectors between neighboring points. The neighbors are on the boundary between successive row sweeps.
dex, dey
are likewise vectors between neighbors of successive zone pairs during one sweep.
p0, p1
are pressures which are calculated as the reciprocal of the zone volume, Boyle's law, with a factor of dt thrown in to save a few subsequent multiplies.
ddx, ddy
is the vector along the boundary between the two triangles of a cell.

Dynamics

Every time cycle each vertex is pushed on by the six triangles that meet that vertex. The push is proportional to the pressure and the side (as vector) of the triangle opposite the vertex. The vector is rotated by 90°. The push is manifest as an increment to the momentum or velocity of the vertex. In this code the mass of a vertex is taken as 1.

Each cycle the position of the vertex is incremented according to its velocity.

Initial Configuration

Here we use the notation (i,j) to denote the vertex described by a[i][j], There are two populations of triangles. For 0≤i<m and 0≤j<n vertices, (i,j), (i,j+1), (i+1,j) describe a triangle of the first population and (i,j+1), (i+1,j), (i+1,j+1) describes a triangle of the other population. These are all the triangles and they are all initially equilateral. 
1 2nd subscript varies fastest
2 Subscripts are [i][j]
3 We traverse mesh in cache friendly fashion.

    a[0][1]  a[1][1]
a[0][0]  a[1][0]

a[0][0] is lower left of mesh.

Debugging

Lines 19 thru 29 verify that symmetry of the initial conditions is maintained. This found most of the bugs since the code is much less symmetric than the configuration.

A movie sending binary data thru a channel to a Java program. (incomplete)
Transcribed to OCaml with movies.