# How to implement adaptive mesh refinement using conformal triangles

I am trying to implement adaptive mesh refinement for a finite element code. The code uses (at least for now) linear triangles and so when I do the mesh refinement I want the triangular mesh to remain conformal, i.e. no hanging nodes:

I am planning to use the longest-edge bisection method by Rivara [33,34]. Unfortunately I have not been able to access these papers since I have not found a freely downloadable pdf. I have looked at a couple other papers/resources that describe this method, i.e. fea8.pdf as well as Wolfgang Bangerth's video series lecture 14 - 18, however I am having trouble with the details on how this algorithm would actually be implemented.

What I have done so far is assume an initial coarse mesh, i.e. something like:

In my program my mesh is represented by an array of binary trees where each original coarse triangle is the root node of one of the binary trees in the array. So for example in the above mesh I would have an array of binary trees of size two. Graphically my trees would look like:

Now say for example triangle $e1$ is flagged for refinement. We bisect this triangle by the longest edge to create two new sub triangles:

and tree structure:

The problem is that triangle $f1$ has a hanging node, and worst still the root node in the tree for triangle $f1$ doesn't "know" it has a hanging node along that edge nor does it "know" that this hanging node has number 5. One solution that I came up with would be to have each node representing a triangle to also have additional pointers that point to the neighbours of that triangle. This however quickly makes the array of trees really complicated with pointers between nodes all over the place trying to keep track of which triangles are neighbours of other triangles.

In C++ for example my binary tree data structure looks something like:

template <class T>
class Tree
{
private:
struct Node
{
T element;
Node *parent;
Node *left;
Node *right;

Node *neigh01;
Node *neigh12;
Node *neigh20;

int hnode01;   //store hanging node numbers for this node
int hnode12;
int hnode20;

bool refine;   //refinement flag saying this node is in need of refinement

Node(T e) : element(e), parent(NULL), left(NULL), right(NULL), neigh01(NULL), neigh12(NULL), neigh20(NULL), hnode01(-1), hnode12(-1), hnode20(-1)
{
}
};

//rest of class ...constuctors, methods, etc...
};


where you can see I have included additional node pointers that attempt to keep track of where neighbor nodes are stored. Actually implementing this bisection method for adaptive mesh refinement appears to be quite complicated and I am not sure how to go about doing it. Is my idea of having additional pointers to keep track of neighbours the right way to go about solving this? Does anyone have freely available sources that explain the bisection algorithm in more detail, i.e. explaining the data structures used etc? Can anyone point me to an algorithm for adaptive mesh refinement using bisection? The one written in fea8.pdf seems to indicate that we first run through all nodes that need refinement and split them. We then go through (iteratively?) to correct for hanging nodes. Is this correct/what is usually done? Is my idea of represent my mesh as an array of binary trees a good/standard way of representing a mesh that will be adaptively refined? Basically I need any help I can get in implementing this adaptive mesh refinement algorithm.

Thanks

[33] M.C. Rivara. Design and data structures of a fully adaptive multigrid finite element software. ACM Transactions on Mathematical Software, 10:242 264, 1984.

[34] M.C. Rivara. Mesh refinement processes based on the generalized bisection of simplices. SIAM Journal on Numerical Analysis, 21:604 613, 1984.

• Note that I don't write this sort of code professionally or anything, so I can't say for sure if you're on point. However, I think the tree structure makes sense for the triangle. Have you considered doing a similar thing for the edges? Then a given triangle could check it's edges and make sure they still have null for the children. If the edges had children, then you would know there's a hanging node. – spektr Oct 11 '15 at 21:19
• I haven' programmed any meshing except for simple quadrilaterals, so I'm not that knowledgable. Why not create a connection matrix which will be sparsely populated and can be boolean. At this level the problem gets similar to some graph analysis methods used for solving electrical networks. You can create it so it has nodes$\times$nodes, or nodes $\times$triangles as axes. Another possibility is to add a dimension, or put in pointers to your binary trees? In the refinement step you would expand those matrices and fill it further up. I hope it helps in some way! – WalyKu Oct 15 '15 at 8:02
• @choward After thinking on this problem further I think that I will in fact need each node in the tree representing a triangle to have additional node pointers. These pointers will point to other neighboring nodes in possibly completely different binary trees. As each node is refined by adding a left and right child I must update/assign the node pointers in the child nodes to point to the child nodes' neighbors (which might exists in separate binary trees). I think a similar thing would have to be done with your edge idea as each edge needs to know what triangles it is a part of. – James Oct 18 '15 at 6:34
• Well, as of Monday, I came into need of exactly the same thing, but this is the first time I am doing something like this; although it does not look difficult, I could use a bit more details of the algorithm: For example, in the example above, if point 2 was lower twice as far from point 3, how is this solved? I mean, say cell 134 was divided first and point 5 now exists; then (with point 2 further down) segment 1-3 is not the longest edge of cell 123...so? Is cell 123 divided within a single visit until it uses point 5? At a more global level, is the mesh re-built after every single pass? and – user21974 Oct 12 '16 at 17:59
• @gsal Yes if node 2 is twice as far down from node 3 then things do become more complicated. What you must do is run through the mesh refining (along their longest edge) the triangles that are marked. Doing this may create situations as you describe where there are hanging nodes. You must then iteratively go through the mesh again (refining as needed) to remove all hanging nodes. – James Oct 14 '16 at 5:27

The way this is typically handled is by running a pass before you actually split cells that determines which other cells also have to be refined. In pseudo-code, this could look like this:

for (cell in cells)
if (cell->is_marked_for_refinement && cell->edge_to_split == unassigned)
{
e = longest edge of cell (0...2);
cell->set_edge_to_split(e);
n = cell->neighbor(e);
n->mark_for_refinement();  // if we split this cell along edge e,
// then we also need to split cell n
n_of_n = n->get_neighbor_index(cell); // find out along which edge
// of n the current cell is
n->set_edge_to_split (n_of_n);
}


This process ensures that you first mark all cells that actually need to be split, before you do the splitting.

• I think my problem is how to determine the neighbor cells efficiently, i.e. in the above pseudo-code the line: n = cell->neighbor(e);.If I understand your answer correctly then each cell (or maybe each edge) has pointers to its neighbor cells, but as the mesh is refined these pointers need to be updated to point to the smaller cells being created. This needs to be done in some efficient way and I am not sure how to do this. – James Oct 12 '15 at 17:13
• Correct. You will need neighbor pointers for many other operations as well. You can efficiently set the neighbor pointers of new cells if you have neighbor pointers between their parents. After you have created all new cells, you essentially need a loop over all (now no longer active) parent cells and then do something like parent->child(c)->set_neighbor(n, parent->neighbor(n)->child(c) for some permutation of indices c=0,1 and n=0,1,2. – Wolfgang Bangerth Oct 12 '15 at 19:22

With the help of Wolfgangs answer I think I have figured it out. It is necessary to include additional pointers in my Node class so that each Node in the binary tree can know who are its neighbors. For example my Node and Tree classes look like:

template <class T>
class Node
{
public:
T element;
Node *parent;
Node *left;
Node *right;

private:
Node *neigh01;
Node *neigh12;
Node *neigh20;

int hnode01;   //store hanging node numbers for this node
int hnode12;
int hnode20;

bool refine;   //refinement flag saying this node is in need of refinement

public:
Node(T e) : element(e), parent(NULL), left(NULL), right(NULL), neigh01(NULL), neigh12(NULL), neigh20(NULL), hnode01(-1), hnode12(-1), hnode20(-1)
{
}

//rest of methods
};

template <class T>
class Tree
{
private:
int numNodes;
int numLeafNodes
Node<T> *root;
Node<T> *iterator;

public:
Tree();
Tree(T e);
~Tree();

void insert(T e);
void deleteTree(Node<T> *leaf);

//rest of class methods, etc...
};


The three additional Node pointers - neigh01, neigh12, & neigh20 - allow us to keep track of any neighbor Nodes. Graphically our tree structure might look like:

where the red lines indicate neighbor Node connections. The important thing to remember when a triangle is refined (i.e. split by longest edge bisection) is that these new child Nodes can only be neighbors with:

1. NULL (in the case of a boundary, i.e. no neighbor)
2. each other (each left child Node will always have a neighboring right child Node)
3. neighbors or children of the child Nodes' parent Node (any Nodes that are neighbors or children of a parent Node may be neighbors of the newly created child Nodes)

It is a little tricky, but with some shallow recursive functions I was able to ensure that the correct neighbor pointers were assigned every time a Node is refined. Here is the result for a test mesh before and after refinement:

Indeed the mesh remains conformal!

• I have been trying to implement a red-green adaptive refinement algorithm but couldn't reach the right implementation. Could you please explain your longest edge-based refinement algorithm implementation that would give me an idea to implement a red-green refinement algorithm. Thank you. – Aravindh SK Jan 18 at 15:44