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?


1 Answer 1


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.

  • $\begingroup$ Thanks Wolfgang for your answer, and also for the last paragraph. Do you have some references about the technique you highlighted in the first two paragraphs? I don't want to reinvent the wheel, but I was trying to implement P2 basis with a 2D domain since I'm actually using dealii in my course, but the assembly part always triggers me @WolfgangBangerth $\endgroup$
    – FEGirl
    Mar 5, 2021 at 21:53
  • $\begingroup$ 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 $\endgroup$
    – FEGirl
    Mar 5, 2021 at 22:45
  • $\begingroup$ 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. $\endgroup$ Mar 6, 2021 at 1:38
  • $\begingroup$ The enumeration of nodes and vertices is conceptually entirely separate. On every cell you could have an array of length 6 that stores the indices of the degrees of freedom that live on that cell. You have some sort of convention where these are located -- for example, the first 3 are on the vertices of the cell, the next 3 are on the edge midpoints. Then you just need to make sure that the 6 indices you store on one cell are consistent with the 6 indices you store on the neighboring cell. In essence I just explained to you what DoFHandler does in deal.II :-) $\endgroup$ Mar 6, 2021 at 1:40
  • $\begingroup$ Thanks @WolfgangBangerth, I'm taking a closer look at this design choice $\endgroup$
    – FEGirl
    Mar 7, 2021 at 10:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.