I use a variant of option 3, that is: N independent containers, one for each cell/element type/order such that all cells of the same type are aligned in memory. When I want to perform operation x on all cells (integrate, quadrature, evolve, numerical flux...) I do:
for each container<cell_type> in cell_containers
for each cell in container<cell_type>
do x<cell_type>(cell)
where x is a different, specialized and optimized function for each cell type. To allow random access a small vector of pointers to all cell containers is stored, inside each container random access iterators are available (which should be at least as fast as the second approach but here your jumps in memory are constant inside each container because the elements of the same type are aligned sequentially). If you need global ordering for say a space-filling curve I can't see how implementing it to iterate over elements[order][id] should be more difficult than elements[id] but I haven't tried it.
The main advantage is that the compiler knows a lot about what is going on and can safely optimize (no virtual dispatch, no aliasing, no breaks out of the loop, the types and element order inside your numerical kernels are all compile-time constants...). As a nice side effect, hybridizing an MPI code using this approach is just as easy as loading the std parallel library(*) and compiling with OpenMP, no extra code required. For boundary elements one has usually to perform other operations and store different data. 904
I store them separately using this same approach and define specialized kernels for them.
The main drawback is that you need to perform memory allocation very carefully. The problem is, before hand the only thing you know is how much memory is available, but you don't know how many elements of which order you are going to need. In particular, imagine that you have allocated a certain amount of memory for a certain type of cell, and you want to exceed it. This will require that you have at least enough free memory to copy the whole vector of cells plus one extra cell. It also involves copying all the cells from that vector to the new one but that is fast. This is not smart, you probably want to allocate memory for more than one extra cell, but for how many? You probably dont want to allocate memory for twice as many cells. If you want to do this in a smart perfomant way you really need your own allocator that knows about your data structures, your system, and the on going simulation.
Pros:
* Operations that need to be performed on all cells (numerical kernels) are very fast.
* Random access is available and is as fast as with the "one vector of cell data/one vector of cell ptrs" approach.
Cons:
* Memory allocation has to be performed very carefully (but isn't this always the case?).
* using TBB requires replacing your for_each loops with parallel_for_each or you can just implement your own for_each that takes a parallelization policy (or macro it: please dont do this). 77