Two-dimensional mesh in fem: generating P1, P2, P3,... mesh from a P1 mesh

Question

I have a two-dimensional mesh generated by triangle (the mesh generator software is not relevant). This software generates a perfect mesh for approximate the solution by piecewice linear functions (P1 mesh).

For example, for the unit square [0,1]x[0,1] in 2D I have a file with the coordinates of its nodes (for example, a mesh with 5 nodes):

1   0.0 0.0  # coordinates of node 1
2   1.0 0.0  # coordinates of node 2
3   1.0 1.0  # coordinates of node 3
4   0.0 1.0  # coordinates of node 4
5   0.5 0.5  # coordinates of node 5

called coordinate.dat, which have 4 elements (triangles) whose conectivity is given by the file called element.dat:

1   1 5 4  # vertices of triangle 1
2   1 2 5  # vertices of triangle 2
3   2 3 5  # vertices of triangle 3
4   5 2 4  # vertices of triangle 4

Now, I need to program problems where the finite element spaces belong to arbitrary polynomial degree, for example, approximate the solution by polynomials of degree 2, 3, etc.

To do this, I need to build a new mesh with new nodes and a new conectivity file with more nodes. An example of the general mesh that I need:

I already know how to calculate the coordinate of all new nodes on each element (it is just a convex combination easy to calculate) but I can't get a easy way to generate the new conectivity file.

I've been thinking for several days how to program a general method for this, but I could not come up with any.

Do you know some method to create meshes P2, P3, ... from a mesh P1? or more precisely, How can I build the new conectivity file?

Answer 1

The mesh for Pk elements is usually exactly the same as for P1. The only difference is that you also need to keep track of your edge indices for higher orders and associate degrees of freedom (DOFs) with edges (from P2 on) and elements (from P3 on). This is only a matter of defining local indices for nodes and edges with respect to your reference element, mapping to the actual mesh element and doing some index bookkeeping. Computing indices is quite straightforward if you figured that out and does not involve working with refined mesh data structures. The advantage of this technique is that you do not rely on a specific "placement" of DOFs on mesh nodes and generalizes easily to more general finite elements. I suggest you start over from scratch and consult the existing literature on the implementation of finite element methods. There is even a book with the title "Understanding and Implementing the Finite Element Method" by Gockenbach. While I have not read it myself, I am pretty sure that he explains these basics better than I ever could in this post. Other popular introductions usually also have a chapter on the implementation.

What you intended to do is however not a useless exercise. A popular approach for visualizing higher order polynomials is to use refined meshes, see also this question. In this case you need to deal with an interpolation to a refined mesh. Interpolation is easy if you know the vertices and refining is usually done recursively, i.e., you refine once, then again, etc. In this way it is easy to keep track of connectivities if you have a suitable index for your edges. Setting that up is a straightforward task since edges can be defined as a tuple of node indices which you can (again) define once and for all for your reference element and map to each of the actual mesh elements. However, since your question is related to setting up higher order elements, I will not elaborate on mesh refinement here.