# Programming Finite Element Methods in C++

I am trying to develop a library for finite element methods in C++ and for that I am looking at the data structures for meshes. Based on what I've read up on fenics and deal.ii, the general understanding seems to be that a basic mesh class should have the following members:

1. Cells, Faces, vertices.
2. Connectivity relationships (What are the neighbors of a given cell, common faces between cells, vertices of a given cell, boundary indicators for cells, faces and vertices)
3. Cell iterators

These are my questions:

1. What is the best way to store cells of a mesh. I think an STL vector would suffice. Is this the popular option? I have also read that linked lists have been used for cells but I cannot see the advantage of linked lists over vectors in this case.
2. Again, for the connectivity relationships, STL maps seem to be the best method for me. Any comments on this?
3. Does the same data structure work for both triangular and quadrilateral meshes. I think it should but I am not sure if I am missing anything basic here?
4. Are there any more connectivity relationships to be kept in mind?
5. Are there any good references for programming a finite element library?

Thank you.

• I guess all of your questions may be answered checking source code for the libraries you mentioned. What features are missing in dealii or fenics? I personally have experience only on dealii and they have quite an active community, you really want to work all alone? Oct 29, 2014 at 9:52
• Thank you Nicola. I am working on this project out of interest in building software libraries. I am assuming that dealii and fenics use advanced/efficient programming concepts for the questions I mentioned above and I wanted to know if vectors and maps are good enough for basic finite element libraries.
– gk1
Oct 29, 2014 at 10:02

One of the authors of fenics, A. Logg, have written a very good paper on datastructures of storing meshes. The paper is A. Logg (2009). Efficient Representation of Computational Meshes http://arxiv.org/abs/1205.3081

In fact it's always a tradeoff between storing all the topological informations (nodes around nodes, faces around nodes, etc...) OR having to recompute them before using them. If you store them then you will probably end up with a very massive amount of data to store, if you recompute them using iterators will take some (too much) time. Depending upon the algorithm (typically for FEM, with linear elements you just want nodes of polyhedrons to assemble your local stiffness matrices), you may want to use some specific traversals and then you can generate what you need on the fly.

Anyway, the proposed solution in fenics is based on many arrays of integer handles and is probably what's closest to the most efficient method (both in terms of memory and access).

In deal.II, we basically only use vectors. Maps are too slow and scatter data all around memory, so we typically don't use them if the keys are integers and within a given range. For example, for the connectivity between cells, you can do arrays (STL vectors) in which you store neighbor indices and so that neighbor indices $4i\ldots 4i+3$ correspond to cell $i$ (in 2d, because there are 4 neighbors per cell). This way you can do direct array lookup, rather than walking through the tree in STL maps.

You probably realize this already, but we have put a lot of research into the data structures we use in this library, and the scheme we use has led to exceptionally low numbers of cache misses. You will spend a lot of time coming up with the same kinds of solutions we use if you want to be efficient.

• Would you recommend vector of int or vector of data structures? i.e. std::vector<int> vs std::vector<Vertex&>. Is this relevant to cache misses?
– user1800
Feb 21, 2016 at 3:54
• @JayeshBadwaik, your question is whether it is better to store a structure of arrays vs an array of structures. Our experience is that the former is the better option because in any given algorithm, you typically only use a small subset of the data stored for every cell, but you will walk from cell to cell. Consequently, you want to go through the memory that corresponds to arrays in a linear fashion, rather than selectively reading a few pieces and data and then skipping forward by a large offset. The latter creates cash misses. Feb 21, 2016 at 4:39