How to implement adaptive mesh refinement using conformal triangles

Question

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.

Answer 1

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.

Answer 2

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!