# Why should one use a tree structure to represent discrete function spaces?

In FEM/FV codebases, I stumbled upon the fact that the discretized functionspaces are represented within the code as a tree structure. I find this very puzzling.

Example: lets say somebody wants to simulate the navier stokes equations. The codebase expects the discretized function spaces to be put into a data structure like:

So an hierarchical tree structure. I do not understand why that should be a reasonable choice because:

1. There is no intrinsic hierarchy of the fields. In the image, the pressure variables are listed at a higher level than the velocity components. But that is purely cosmetic display. There is no intrinsic hierarchy of the components influencing our coupled system. If there is no hierarchy, then a tree structure does not seem sensible to me.

2. Data locality. I have an option where I can tell the programm to align the data cell-wise or in lexicographic ordering. That is a nice feature, but it is, at it's heart, the age-old question of (array-of-structs vs. struct-of-arrays). This also does not explain to me why I should use a tree.

3. Limited Complexity. I have not come across any codebase or scenario, where the function spaces are so many, or so complex, that that would justify a non-trivial data structure. Even if you do full MHD simulations, you could easily fit your fields into (array-of-structs vs. struct-of-arrays). And it is questionable whether a tree structure would simplify or confuse your code.

What am I missing? Why would anybody use a tree structure for this? (Thanks in advance!)

• OK... If you think using tree structure doesn't make sense, so what do you think might be better option here? Oct 23 '19 at 16:14
• What you need in your assembly procedures is all the local degrees of freedom. If I were to implement it (as a newby), I would just use a container for all dof's associated with a cell/element. Oct 24 '19 at 7:38

## 1 Answer

Although I can't speak for NS or MHD, I do find this "componentization" of function spaces to be a useful design principle in CEM (computational electromagnetics), especially for high-order (in p) discretizations.

CEM often uses multiple spaces at the same time: grad-conforming functions to represent electric potential, curl-conforming functions to represent E-field, div-conforming functions to represent B-flux and currents, etc. These spaces are related to one another via the deRham diagram, the range of each differential operator (grad/curl/div) forms the kernel of the next (the exact sequence property).

In my experience this can (should?) show up all the way at the implementation level, because it saves you time and helps guarantee correctness. For example, if you've spent time and money tabulating/implementing a high order grad-conforming space (basis functions $$\phi$$'s and their gradients $$\nabla \phi$$'s, you can (should?) reuse those $$\nabla \phi$$'s as basis functions within the curl-conforming space, then you just need to enrich them by adding whatever additional functions are that are needed span the range of curl to your desired order. Similar constructions hold between curl-conforming and div-conforming spaces (the curls of your curl conforming functions should be members of your div-conforming basis).

Unfortunately, this enrichment process is messy/broken down at p=0 because of unisolvency issues (p=0 curls already spans p=0 gradients, etc). It really only shines at p >= 1. In practice, I use the following decomposition:

• hgrad0: p=0 gradient-conforming functions
• hgradp: p>0 gradient-conforming functions
• hcurl0: p=0 curl-conforming functions
• hcurlp: functions that "enrich" the range of curl for p>0
• hdiv0: p=0 div-conforming functions
• hdivp: functions that "enrich" the range of div for p>0

With these building blocks, you can tabulate the whole diagram:

• hgrad := hgrad0 + hgradp
• hcurl := hcurl0 + grad(hgradp) + hcurlp
• hdiv := hdiv0 + curl(hcurlp) + hdivp
• L2 := piecewise constant + div(hdivp)

And each of these spaces (i) splits apart the kernel from the range of the respective differential operator and (ii) splits apart the p=0 space from the p>0 .. both of these properties are useful when designing preconditioners.

• Thank you very much for this answer. I need to do some reading I guess:-) Oct 24 '19 at 7:53
• Lets say the global function space is the following: B_th = B_v1 x B_v2 x _Bv3 x B_p. Then the grouping of them still seems a bit arbitrary. I guess the whole thing comes from the fact that the involved operators only apply to subspaces of the whole problem, and then the grouping is induced! Oct 24 '19 at 8:55