Renumbering the nodes for quadratic basis functions for a 2D domain

Question

I have a simple triangulation for a 2D domain, described by the connectivity matrix $T$ and by the point matrix $P$.
For didactic purposes, I assembled the stiffness matrix for $-Delta u = f$ by using the reference triangle. With $mathbb{P}1$ basis functions everything is easy because all the gradients are constant there.
Now I'm trying to use $mathbb{P}2$ basis function. In this case we have 6 nodes for the reference element: 3 for the vertices and 3 at the midpoints. indeed the local stiffness matrix must be a 6 by 6, and here be dragons. I don't know how to change the matrices $T$ and $P$ in order to describe the mesh.
Do you have any good reference with also some codes? Or, even better, a Python snippet that given the two matrices $P,T$ returns the new matrices considering also the midpoints?

Wolfgang Bangerth · Accepted Answer

A good abstraction would realize that there is a difference between the vertices of the mesh on the one hand, and the nodes of the finite elements on the other. If you conceptually separate these, then the element you are using has no effect on your matrices $T$ and $P$. Rather, you'd have to have additional arrays that indicate which node indices are live on each of the cells, and where on these cells are located. You wouldn't have to store their locations, though, since that can be computed.
For example, if you had an integer array of length 6 on each cell, you would adopt the convention that the first three indices correspond to the nodes located on the three vertices of the cell (in this order, and you can look up these vertices' indices via the $T$ array and then the location via the $P$ array), and the remaining three are the nodes located on the edges of the cell (for which you can get the locations based on the locations of the vertices on this cell).
This is an example of a more general design principle: Things that are conceptually different should be stored separately. The fact that in some cases one can get away with less really just means that if you want to implement a more general case, you have to work so much harder.

Renumbering the nodes for quadratic basis functions for a 2D domain

One Answer

Add your own answers!

Ask a Question