Xiss,

You are right, that coordinate_matrix is very slow for this example. Here is the reason: the reserve() function gets triggered several times (at least 9 with your example code). reserve() itself calls sort(). And every time it has more and more elements (starting with 15000 on the first call, and finishing with about 2.7 million elements on the last call). These calls to sort() become more expensive as the number of elements grows, and so your main loop spends 80% of its time doing the sort() (I actually profiled it!).
I don't think this can be remedied, because you do have lots of elements (numTetrahedrons*12*12 = 2.8 million) to append to the coordinate_matrix. So I guess coordinate_matrix is not appropriate for this use case after all.

BTW, I think you have a bug in the attached code: r*numNodes will most likely overflow and will become smaller than mtrand.max(), which will make ni be equal to 0, so all elements of the index matrix will be 0, 1 or 2. I'm pretty sure that's not what you want. One easy fix is to declare r as double (instead of int).

-Vardan