Quake solver code

MATRIX3 K;

/* matrix and vector assembly */

FORELEM(i) {

...

}

/* time integration loop */

for (iter = 1; iter <= timesteps; iter++) {

MV3PRODUCT(K, disp[dispt], disp[disptplus]);

disp[disptplus] *= - IP.dt * IP.dt;

disp[disptplus] += 2.0 * M * disp[dispt] -

(M - IP.dt / 2.0 * C) * disp[disptminus] - ...);

disp[disptplus] = disp[disptplus] / (M + IP.dt / 2.0 * C);

i = disptminus;

disptminus = dispt;

dispt = disptplus;

disptplus = i;

}