I've asked this question on Computational Science SE too.

Consider the following triangulation of a rectangle $(0,a)\times(0,b)$:

enter image description here

We have two types of triangles:

enter image description here

I enumerate the rectangles(consisting of two triangles) starting at the bottom left corner from left to right and bottom to top. The triangles are enumerated in the following way: Each triangle of type A has the same index as its corresponding cell and each index of a triangle of type B is the sum of the index of its corresponding cell "$+$ the total number of cells".

For example, in the triangulation depicted above, the cell at the bottom left corner has index $1$ and the contained triangles of type A and B have the index 1 and 17, respectively.

The local numbering of the vertices is shown the second figure. Starting from the bottom left corner, from left to right and bottom to top, I build the global indexing of the vertices. For example, in the triangulation depicted above, the vertex at the bottom left corner has index $1$, its horizontal neighbor has index $2$ and its vertical neighbor has index $6$.

This indexing is suitable if each triangle is considered to be a linear Lagrange element (evaluation at the vertices). How do I need to modify my indexing, if I want to use quadratic Lagrange elements instead (evaluation at the vertices plus evaluation at the midpoints of the edges)?

I would like to use the following local numbering for the vertices:

enter image description here

It's not clear to me how I should build the global numbering of the additional nodes (i.e. the midpoints of the edges). I'm not sure, but I could imagine that the choice of the numbering influences the linear system of equations in a crucial way (I've read, for example, that a tuple of increasing local vertex numbers should correspond to a tuple of global vertex numbers).

  • 3
    $\begingroup$ You did not say what sparse solver you are using. All modern sparse solvers have a step where the equations are internally renumbered to minimize fill-in during factorization and so lower the computational cost. In that case, you can choose the global node numbers using whatever scheme you find most meaningful to you. $\endgroup$ Mar 24, 2017 at 16:57
  • $\begingroup$ @BillGreene I don't use any solver. I'm trying to understand a simple implementation for educational purposes. $\endgroup$
    – 0xbadf00d
    Mar 24, 2017 at 19:16
  • $\begingroup$ One simple example, would be Cuthill-McKee algorithm that permutes (renumbers) a sparse matrix for bandwidth minimization. That would allow shifting the 'renumbering issue' to the solve step of the algorithm (where it only matters), not to the fill-part. $\endgroup$
    – Anton Menshov
    Mar 24, 2017 at 22:21

1 Answer 1


First approach: using pattern of sparse matrix

You are talking about triangular Lagrange FEs of second order, so I suppose you are familiar with ″ of first order—you know how to assemble system resulting from Laplace or Poisson equation using linear elements.

For the resulting matrix we typically use compressed sparse (lower triangular) column format, CSlC (see SPARSKIT documentation or User’s Guide for the Harwell–Boeing Sparse Matrix Collection for details).

First we build the pattern of our matrix (vectors colptr and rowind), and then we assemble our system. When one has matrix pattern, one also has numeration for middle nodes (or ribs)!

Consider the following mesh and resulting matrix pattern:

enter image description here

We have $$ \text{colptr} = \left\{ 1, 3, 6, 10, 11, 11, 12, 12 \right\}, \\ \text{rowind} = \left\{ 1, 7, 3, 4, 7, 4, 5, 6, 7, 5, 7 \right\}. $$

In order to get element $(i, j)$ of CSlC matrix we do

for k = colptr[j], colptr[j] + 1, ..., colptr[j+1] - 1
    if (i == rowind[k]) return values[k]

Typically any vertex in triangulation has degree ~6, so this operation is $O(1)$.

So if you want to get a number of a rib connecting vertices $i$ and $j$, simply use $k$ from the algorithm above. In my picture, for one, a rib connecting vertices 3 and 7 has number 9.

Second approach: using elements’ neighbors

You are talking about rather structured mesh. In reality, in order to handle more complex geometries, we typically represent a mesh as vector of vertices + vector of elements (triplets of indices of vertices). Sometimes it is also convenient to store vector of neighbors (triplets of indices of neighboring elements). For one, Mathematica does so in its mesh abstractions (see Scope > "ElementConnectivity"). I bet that all popular software tools for mesh generation such as gmsh can do this as well.

Clearly, if you have neighbors, it is easy to enumerate ribs iterating over elements. I tend to use this approach in my projects since I have access to neighbors.


but I could imagine that the choice of the numbering influences the linear system of equations in a crucial way

Indeed, different numeration of DOFs corresponds to different pattern of matrix. However, one should not care about this at all if one wants to use iterative solver.

In contrast, if one wants to use some kind of sparse direct solver based on LU decomposition, one should care about DOFs numeration since it affects bandwidth (hence, memory) crucially: LU decomposition of CSC-type matrix may be dense; LU decomposition of Skyline matrix preserves pattern of sparsity.

However, as @BillGreene and @anton-menshov mentioned, there are algorithms that minimize bandwidth of your matrix (i.e. renumber DOFs). You should not care about this while assembling / choosing DOFs numeration.


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.