Creating a custom high-order Lagrangian element with Hierarchical shape functions using basix

I read the documentation for Basix and it seems that Lagrangian elements using Hierarchical shape functions are not included among the pre-provided element types. So I decided to try creating a custom element.

But it seems there is no way to attach multiple DOFs to one node. Take the triangular element based on Hierarchical shape functions as an example:

As shown in the figure, using Hierarchical shape functions can increase the element order without changing the original shape functions.
But I refer to the example in Creating a custom element — Basix 0.9.0.0 documentation which seems unable to create such the element. I’ve tried many combinations, adding elements of wcoeff and x, changing the size of M, and all produce dimension errors, or NaN output. If I’m missing something please let me know, I’d be very grateful.

To make it clear what form of shape function interpolation I want to implement, I attached my code. Obviously this code contains errors.

import basix
import numpy as np
from basix import CellType, MapType, PolynomialType, LatticeType, SobolevSpace

wcoeffs = np.zeros((9, 6))

wcoeffs[0, 0] = 1
wcoeffs[0, 1] = -1
wcoeffs[0, 2] = -1

wcoeffs[1, 2] = 1
wcoeffs[1, 4] = -1
wcoeffs[1, 5] = -1

wcoeffs[2, 1] = 1
wcoeffs[2, 3] = -1
wcoeffs[2, 4] = -1

wcoeffs[3, 2] = 1

wcoeffs[4, 2] = -1
wcoeffs[4, 5] = 1

wcoeffs[5, 4] = 1

wcoeffs[6, 1] = 1

wcoeffs[7, 1] = -1
wcoeffs[7, 3] = 1

x = [[], [], [], []]
x[0].append(np.array([[0.0, 0.0]]))
x[0].append(np.array([[0.0, 0.0]]))
x[0].append(np.array([[0.0, 0.0]]))

x[0].append(np.array([[1.0, 0.0]]))
x[0].append(np.array([[1.0, 0.0]]))
x[0].append(np.array([[1.0, 0.0]]))

x[0].append(np.array([[0.0, 1.0]]))
x[0].append(np.array([[0.0, 1.0]]))
x[0].append(np.array([[0.0, 1.0]]))

x[2].append(np.array(np.zeros((0, 2))))

for _ in range(3):
    x[1].append(np.zeros((0, 2)))


M = [[], [], [], []]

M[0].append(np.array([[[[1.]]]]))
M[0].append(np.array([[[[0.]]]]))
M[0].append(np.array([[[[0.]]]]))

M[0].append(np.array([[[[1.]]]]))
M[0].append(np.array([[[[0.]]]]))
M[0].append(np.array([[[[0.]]]]))

M[0].append(np.array([[[[1.]]]]))
M[0].append(np.array([[[[0.]]]]))
M[0].append(np.array([[[[0.]]]]))

M[2].append(np.zeros((0, 1, 0, 1)))

for _ in range(3):
    M[1].append(np.zeros((0, 1, 0, 1)))

element = basix.create_custom_element(
    CellType.triangle, [], wcoeffs, x, M, 0, MapType.identity, SobolevSpace.H1, False, 1, 2)

points = np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
print(element.tabulate(0, points))

I would like finding the right way to create this type of higher order element.