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

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?

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).
• I'm trying to the example on your second paragraph, but I'm stuck here: 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). In order to fill those new three components of $T$ I need to give a numbering to those nodes on the edges, and to me it's not clear how I should enumerate them Mar 5 '21 at 22:45
• I don't have a reference (just 25 years of experience :-) ) but if you want something to look at just look at the separation of the Triangulation and DoFHandler classes: The former essentially stores what you denote by the $P$ and $T$ matrices. The latter deals with degrees of freedom defined on the former. Mar 6 '21 at 1:38