3
$\begingroup$

I have a 2D mesh of triangles used in Finite Element method to discretize the domain. I want to calculate the midpoints of all the edges because I want to use $\mathbb{P}^2$ elements. I am using MatLab. The data that I have is:

  • Matrix containing nodes coordinates, $x$ and $y$. The size is $N\times 2$, where $N$ is the number of nodes.

  • Matrix containing each triangle nodes. The size is $NT\times 3$, where $NT$ is the number of triangles in the mesh and 3 columns because each triangle has three nodes. For example, a row of this matrix could have $[10, 4 ,22]$ so the triangle is formed by these nodes and in order to know their coordinates I just take the rows $10, 4$ and $22$ from the node coordinate matrix (the previous one).

  • Array with boundary nodes.

I know in advance that the number of nodes (including midpoints) is going to be $2N\times N_b -3$ where $N_b$ is the number of nodes in the boundary of the mesh, so I can preallocate the matrix. The new node matrix will contain the previous nodes coordinates and the new ones and the element matrix will now have 6 columns because one node will be added per edge.

What my code is doing is to run across all triangles and calculate their edges midpoints and if it is not yet in the coordinate node matrix include it (here I am losing many time checking if it is yet inside or not). The main problem is that as interior edges are common to two triangles many times the midpoint of an edge is already calculated so I am wasting time calculating the midpoint and after it checking that is already in the matrix.

Could you give me any advice about how to calculate more efficiently the midpoints of edges and include them in the matrix. Note: I do not have a matrix with the edges of the mesh.

Thanks!

Improve this question
$\endgroup$
2
  • $\begingroup$ As a possible alternate solution. Use gmsh to create your mesh and then just write a small script to import the mesh into matlab. Its not too hard as I have done this myself. This will allow alot more flexibility i think in terms of different kinds of elements and order of elements. $\endgroup$
    – James
    Jun 22, 2015 at 16:03
  • $\begingroup$ In this case I am using Distmesh but it only returns edges at the boundary. The main advantage of Distmesh is that I can create the mesh automatically from my code. However, Gmsh does not allow it easily. But thanks for the suggestion. $\endgroup$
    – James
    Jun 22, 2015 at 18:05

2 Answers 2

Reset to default
4
$\begingroup$

If you have an array that stores the indices of the 3 neighbors of each cell, then you would only compute the midpoint of an edge of the neighbor cell has a higher index than the current cell, or if there is no neighbor at all. This way you have an easy tie breaker to decide which of the two cells is responsible for computing the edge midpoint.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I understand that you are suggesting to use like a 'protocol' so no edge midpoint is calculated twice. However, I do not get the idea of the 'protocol'. Do I need to know the index of the neighbor cell? Where is this information taken from? $\endgroup$
    – James
    Jun 22, 2015 at 18:13
  • $\begingroup$ The 'index' of the cell is the number of the row of your $N_T\times 3$ matrix in which you store the indices of the three vertices of the $N_T$ cells. $\endgroup$
    – Wolfgang Bangerth
    Jun 24, 2015 at 11:43
  • $\begingroup$ Yes, I know, but do I need to know the index of the neighbor cells? $\endgroup$
    – James
    Jun 24, 2015 at 13:48
  • $\begingroup$ Yes. You need to build that map from every cell to its 3 neighbors. You'll need it in many more contexts in the finite element method. $\endgroup$
    – Wolfgang Bangerth
    Jun 24, 2015 at 14:13
2
$\begingroup$

If the main bottleneck is computing a list (well, in MATLAB, I guess it's a matrix) of unique edges, as long as you can find a hash table implementation, you should be able to find a(n average-case) linear time algorithm in the number of vertices (or elements, or edges, by Euler's formula, since your 2D mesh is a planar graph).

On way to do this is to iterate over the matrix containing the triangle nodes. Each row contains three node indices, and thus three two-element vectors of potential edges. Insert these vectors into your hash table, using the sorted version of the vector as a key. (An example: If my edges are $[1, 3]$, $[3, 2]$, $[2, 1]$, my hash table should look like {[1, 2]: [[2, 1]], [1, 3]: [[1, 3]], [2, 3]:[[3, 2]]}, where I'm using Python dictionary notation to represent the hash table because Python dictionaries are implemented as hash tables in the standard library.) Then reading off the keys of the hash table will give you unique edges, and you can use this information to insert midpoints.

There are almost certainly more efficient data structures for manipulating finite element information. Sieves, which I think is DMPlex in PETSc, come to mind; other DAG-like representations exist in the literature that are similar in spirit, if not practice. However, I don't think there's a MATLAB library that implements sieves, whereas there are hash table implementations available for most languages as user-facing libraries, making it an easier solution to kludge together in any language. (It also depends on less in the way of data structures and algorithms knowledge.)

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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