I am trying to construct a hexagonal/geodesic grid on a sphere. For my purposes I am only focussing on the North Pole. So far, I have managed to construct a triangular mesh based on the icosahedron, using the `stripy` package. The `stripy` package allows me to refine the grid through bisection: each edge is split in half, or equivalently, each triangle is split into 4 smaller triangles.

![](https://i.stack.imgur.com/nSWGJ.jpg)
![](https://i.stack.imgur.com/RbEYC.jpg)

I want to create an (almost) hexagonal grid by combining 6 triangles (5 at the pole) as follows:

![](https://i.stack.imgur.com/j1euE.jpg)

I have the following information available:

* latitude/longitude of all vertices as `np.array`
* triplets of indices of triangle vertices

A constraint is that indices are rather 'random', i.e. they don't increase by going outwards in a spiral or something similar.

An option is to find the midpoints of each hexagon(/pentagon) and to group together all the triangles that have this midpoint as one of their three vertices, but I am not sure how to algorithmically go about this. What would be an efficient way of finding the vertices that mark the midpoints of each hexagon? Could it be related in some way to the coarser version of the mesh (i.e. before bisection)?