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!
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:
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.
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)
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.
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.