def METHOD_NAME(
        ldof : ti.int32, 
        NQ : ti.int32, 
        node : ti.template(),  # (NN, 2)
        cell : ti.template(),  # (NC, 3)
        ws : ti.template(), # (NQ, )
        R : ti.template(), # (NQ, ldof, 3)
        S : ti.template() # (NC, ldof, ldof)
        ):
    """
    Note:
        三角形网格上的 p 次拉格朗日元, 组装刚度矩阵
    """
    NC = cell.shape[0]
    NQ = R.shape[0]
    ldof = R.shape[1]
    for c in range(cell.shape[0]):
        x0 = node[cell[c, 0], 0]
        y0 = node[cell[c, 0], 1]
        x1 = node[cell[c, 1], 0]
        y1 = node[cell[c, 1], 1]
        x2 = node[cell[c, 2], 0]
        y2 = node[cell[c, 2], 1]
        l = (x1 - x0)*(y2 - y0) - (y1 - y0)*(x2 - x0) 
        gphi = ti.Matrix([
            [(y1 - y2)/l, (x2 - x1)/l], 
            [(y2 - y0)/l, (x0 - x2)/l],
            [(y0 - y1)/l, (x1 - x0)/l]])
        l *= 0.5
        for i in range(ldof):
            for j in range(ldof):
                S[c, i, j] = 0.0
                for q in range(NQ):
                    vix = R[q, i, 0]*gphi[0, 0] + R[q, i, 1]*gphi[1, 0] + R[q, i, 2]*gphi[2, 0]
                    viy = R[q, i, 1]*gphi[0, 1] + R[q, i, 1]*gphi[1, 1] + R[q, i, 2]*gphi[2, 1]
                    vjx = R[q, j, 0]*gphi[0, 0] + R[q, j, 1]*gphi[1, 0] + R[q, j, 2]*gphi[2, 0]
                    vjy = R[q, j, 1]*gphi[0, 1] + R[q, j, 1]*gphi[1, 1] + R[q, j, 2]*gphi[2, 1]
                    S[c, i, j] += ws[q]*(vix*vjx + viy*vjy)
                S[c, i, j] *= l