The only substantial merit of this code is that it is small and modular. It is otherwise difficult to adapt it to particular jobs.

In general the code is organized to match the ideas in my head. As a result there is a great deal of passing vectors and tensors to and from small routines by value rather than by reference and many of the routines are called at just one place. I am surprised how fast it runs. I have tried to make the math efficient but paid little attention to making the code do the math efficiently. I am tempted to code it for the G4’s AltiVec SIMD. That might approach 3 orders of magnitude speed improvement!

Some unabstracted Data Types defined in Head.h

vec is merely a 3 vector and tri is a 3 by 3 matrix, usually a tensor.

vert represents a vertex. Note that it doesn’t know which zones are its neighbors. Field wh is the location of the vertex and mom is its momentum (= velocity * mass), which is field mass.

lt is a lower triangular matrix:
a00
bc0
def

simp (for “simplex”) is all about a zone. After initialization sh is the zone’s constant shape memory. It is lower triangular from having undergone Gram-Schmidt. This saves a little space and time. The array v holds the indices of the vertices of this zone. mat is the material index, which is constant in this demo.

A link merely denotes a pair of vertices to draw a line between when it is time to draw the configuration.

u is the array of vertices, indexes into which are used by zones to identify their vertices. s is the array of zones. segs is the array of line segments that are drawn to depict the computation.

The Math

This file is nearly purely mathematical.

sm is merely a convenient debugging routine to verify that some value is small.

gs returns the Gram Schmidt ortho-normal vector set.

rx expresses three vectors in the Gram Schmidt coordinates, thus performing a rotation that leaves the matrix in lower triangular form.

shape takes three zone edges in the unstrained geometry and produces the zone’s shape memory as a lower triangular matrix convenient for later processing. This form is the three edge vectors, expressed in an ortho-normal coordinate system that requires only a lower triangular matrix to describe.

cvert creates a vertex but does not register it in the vertex array.

cvec merely packages three floats into a vector. This would not be necessary with the GNU C extensions.

cdvec returns the edge between two vertices expressed as a vector.

csimp together with cxsimp creates a zone (simplex) whose vertices are specified elements from the vertex array. The shape of the zone is computed based on the vertex coordinates, and the zone registered in the array of zones. The mass of the zone is divided equally among its four vertices.

tent registers three zones that form a triangular prism. Caveat: A tent has three quadrilaterals faces. Such a quadrilateral need not be planar. Even if they are initially planar they will not generally remain planar. In fact these faces are divided into two triangles. In order that the zones partition the space it is necessary that two tents that share quadrilateral faces arrange to fold those faces on the same diagonal, thus to share triangular faces, lest some of the space belong to zero or two zones. In tent(i,j,k,l,m,n)

The demo inadvertently conforms to these rules. Be careful—they are not as symmetric as you might hope.

The Graphics

The only thing that I will say about the graphics code (routine draw in Bounce.cp) is that it is short and was written after reading a minimal amount of “Inside Mac”. There seems not to be an easy answer to putting the pixels back upon termination. Evidently I neglect to acquire proper license to write on the screen. The good news is that the code has remained stable since the early Mac SE/30.

Much later (2003): It has now been quite a while since I have tried to run the code. It now shows no signs of running. Here is a Mac Os X version using OpenGL. I could include this code to choose the triangles to be rendered.

main and Initialization

Ad-hoc and mortal code to be thrown out. Even less presentable than the rest. But if you must know the structure zb is the logical form of the vertices. It is an array of slices, each of which locate the vertices of the polygonal slice of the rod. Its elements are shorts (16 bit integers) which serve to index into the array u of vertices where the floating point versions of the dynamic variables live. This plan is predicated on there being a regularity to the zones and vertices.

The Physics

Well there are two overloaded operators defined here in C++ style which are really math but are indeed local to this module.

shove just contributes some momentum to a specified vertex.

push computes the strain tensor εij from the current positions of the three edges of the zone. (Actually it produces εij + δij but δij is deducted just after the call to push).

compress is the compression (actually expansion) undergone by the zone.

The stress is proportional to the strain except for the pressure like component due to compression and the bulk modulus.

The variables frc are computed for each of the three faces adjacent to vertex zero. This is the force due to the stress in the zone upon a zone’s face. It may seem strange to apply that force to the vertex opposite the face but it works out the same in the end and this is faster.
forth is the force for the other face and vertex. The four forces add to zero.

After all of the zones have been processed this way the vertex velocities have been modified and we move the vertices according to their new velocities.

Graphics Code

So far the code alternately computes one δt of physics and paints a screen. With a 400 MHz machine and LCD display this does about 80 frames/sec. The bulk of the time is either in pixel pushing or in mere waiting for the screen to paint. I presume that glutSwapBuffers() sleeps enough to let each buffer get to the screen at least once. Unlike the demo from Apple, there is no display list. Thus deliverFloats must be called for each new frame where the viewpoint has changed or time has passed. The normals are recomputed in deliverFloats even though they do not depend on the viewpoint.