Impulses, FEM and Verlet

After some time struggling with formulations of many diverse categories, I have eventually managed to implement a computer program where the scheme explained in the november entry ( A proposal on a course on Real-Time Structural Dynamics) becomes clear and useful. The main point here is to explain how, by organizing the numerical methods according to what they really discretize, it is possible to grasp and understand such an entangled area of knowledge. One of the main difficulties that I am encountering through this research period is not just that of the thick mathematical language employed but mainly that of mixed concepts when authors come to explain each of the endless methods available in literature. The following is an step by step procedure where time, constraints and material properties are respectively integrated by means of Verlet, Impulse formulation and Finite Elements numerical methods. The software framework for development and visualization is that presented in April 2010 (one year already uuf!): First Ogre+Verlet+Gauss-Seidel simulation. On top of it I have made the necessary modifications and also the adaptation from Linux to Windows so it is easier to make further comparisons with commercial codes. The first that is needed is a set of nodes with their coordinates and their masses, along with the list of constraints that topologically relate one to another and to the environment:

Time integration: Verlet ODE integration method is known to have certain stability but not a big deal of accuracy. However, is fairly easy to implement. An exhaustive explanation on how it works can be found here: http://en.wikipedia.org/wiki/Verlet_integration. In our implementation, it is divided into functions: AccumulateAccelerations and Verlet.
In AccumulateAccelerations we iterate through each mass and add up earthquake readings from a file (per millisecond), gravity acceleration (980 cm/s2), and the internal forces caused by deformation divided by the weight of the mass.
The Verlet procedure, simply updates the position of each particle according to the formula: 

x(t)=x(t-1) + v(t)·dt + a(t)·dt^2

Constraint integration: The only constraints included in the current implementation are those of the distance between two particles. By means of the concept of impulse, once two particles move close or apart by the effect of the accelerations, a corrective force is applied on each one of them. This force is applied in a very short lapse of time, hence can be regarded as an impulse. A more detailed explanation can be found in the work by Jan Bender (www.impulse-based.de). The algebraic system of equations that arises so as to satisfy all of the constraints simultaneously is solved by means of a Gauss-Seidel iterations.

Matter integration: For the consideration of matter properties the program uses the Finite Element Method. A quadratic shape function is used according to the implementation taught in http://www.colorado.edu/engineering/cas/courses.d/IFEM.d/. Chapters 20 and 21 explain clearly how to infer the values for the case of an elastic beam.
This provides our system with a rigidity matrix that is easily attached to each of the distance constraints. As the displacements have already been obtained previously, the solution of the force vector is straightforward by means of a simple rigidity matrix-displacements vector operation:

[f] = [K] · [u]

This force vector stores the effect of the deformation into the constraint and is later used to obtain the particle acceleration in the next timestep.