def METHOD_NAME(self):
    NC = self.mesh.number_of_cells()
    NN = self.mesh.number_of_nodes()
    node = self.mesh.entity('node')
    cell = self.mesh.entity('cell')
    idxi = cell[:, 0]
    idxj = cell[:, 1] 
    idxk = cell[:, 2] 
    v0 = node[idxk] - node[idxj]
    v1 = node[idxi] - node[idxk]
    v2 = node[idxj] - node[idxi]
    area = 0.5*(-v2[:, [0]]*v1[:, [1]] + v2[:, [1]]*v1[:, [0]])
    l2 = np.zeros((NC, 3), dtype=np.float)
    l2[:, 0] = np.sum(v0**2, axis=1)
    l2[:, 1] = np.sum(v1**2, axis=1)
    l2[:, 2] = np.sum(v2**2, axis=1)
    l = np.sqrt(l2)
    p = l.sum(axis=1, keepdims=True)
    q = l.prod(axis=1, keepdims=True)
    mu = p*q/(16*area**2)
    c = mu*(1/(p*l) + 1/l2)
    val = np.concatenate((
        c[:, [1, 2]].sum(axis=1), -c[:, 2], -c[:, 1],
        -c[:, 2], c[:, [0, 2]].sum(axis=1), -c[:, 0],
        -c[:, 1], -c[:, 0], c[:, [0, 1]].sum(axis=1)))
    I = np.concatenate((
        idxi, idxi, idxi,
        idxj, idxj, idxj,
        idxk, idxk, idxk))
    J = np.concatenate((idxi, idxj, idxk))
    J = np.concatenate((J, J, J))
    A = csr_matrix((val, (I, J)), shape=(NN, NN))
    cn = mu/area
    cn.shape = (cn.shape[0],)
    val = np.concatenate((-cn, cn, cn, -cn, -cn, cn))
    I = np.concatenate((idxi, idxi, idxj, idxj, idxk, idxk))
    J = np.concatenate((idxj, idxk, idxi, idxk, idxi, idxj))
    B = csr_matrix((val, (I, J)), shape=(NN, NN))
    return (A, B)