def METHOD_NAME(mesh, surface, gamma=1, theta=0.1): 
    NN = mesh.number_of_nodes()
    NC = mesh.number_of_cells()
    node = mesh.entity('node')
    cell = mesh.entity('cell')
    cellRho, normal= curvaturedensity(mesh, surface, gamma, theta)
    ec = mesh.entity_barycenter('edge')
    cell2edge = mesh.ds.cell_to_edge()
    c = updatecircumcenter(mesh) 
    nodeArea = np.zeros(NN, dtype=np.float)
    massCenter = np.zeros(node.shape, dtype=np.float)
    # 0-th vertex
    a0 = np.cross(ec[cell2edge[:, 2]] - node[cell[:, 0]], c - node[cell[:, 0]])
    a1 = np.cross(c - node[cell[:, 0]], ec[cell2edge[:, 1]] - node[cell[:, 0]])
    a = 0.5*cellRho*(np.sqrt(np.sum(a0**2, axis=1)) + np.sqrt(np.sum(a1**2, axis=1)))
    p = a.reshape(-1, 1)*(node[cell[:, 0]] + c + ec[cell2edge[:, 2]] + ec[cell2edge[:, 1]])/4
    np.add.at(nodeArea, cell[:, 0], a)
    for i in range(3):
        np.add.at(massCenter[:, i], cell[:, 0], p[:, i])
    # 1-th vertex
    a0 = np.cross(ec[cell2edge[:, 0]] - node[cell[:, 1]], c - node[cell[:, 1]])
    a1 = np.cross(c - node[cell[:, 1]], ec[cell2edge[:, 2]] - node[cell[:, 1]])
    a = 0.5*cellRho*(np.sqrt(np.sum(a0**2, axis=1)) + np.sqrt(np.sum(a1**2, axis=1)))
    p = a.reshape(-1, 1)*(node[cell[:, 1]] + c + ec[cell2edge[:, 2]] + ec[cell2edge[:, 0]])/4
    np.add.at(nodeArea, cell[:, 1], a)
    for i in range(3):
        np.add.at(massCenter[:, i], cell[:, 1], p[:, i])
    # 2-th vertex
    a0 = np.cross(ec[cell2edge[:, 1]] - node[cell[:, 2]], c - node[cell[:, 2]])
    a1 = np.cross(c - node[cell[:, 2]], ec[cell2edge[:, 0]] - node[cell[:, 2]])
    a = 0.5*cellRho*(np.sqrt(np.sum(a0**2, axis=1)) + np.sqrt(np.sum(a1**2, axis=1)))
    p = a.reshape(-1, 1)*(node[cell[:, 2]] + c + ec[cell2edge[:, 1]] + ec[cell2edge[:, 0]])/4
    np.add.at(nodeArea, cell[:, 2], a)
    for i in range(3):
        np.add.at(massCenter[:, i], cell[:, 2], p[:, i])
    massCenter /= nodeArea.reshape(-1, 1)
    mv = massCenter - node
    mv -= np.sum(mv*normal, axis=-1).reshape(-1, 1)*normal
    return node + mv