def METHOD_NAME(mesh, surface, gamma=1, theta=0.1):
    NN = mesh.number_of_nodes()
    node = mesh.entity('node')
    cell = mesh.entity('cell')
    normal = surface.unit_normal(node)
    n = np.sum(normal[cell], axis=1) 
    v01 = node[cell[:,1],:] - node[cell[:,0],:]
    v02 = node[cell[:,2],:] - node[cell[:,0],:]
    nv = np.cross(v01,v02)
    length = np.sqrt(np.sum(nv**2, axis=1))
    nv = nv/length.reshape(-1, 1)
    area = length/2
    cellRho = np.abs(9-np.sum(n**2, axis=1))/area
    cellRho = cellRho + theta*np.max(cellRho) 
    sumArea = np.bincount(cell.flat, weights=np.repeat(area, 3), minlength=NN)
    for i in range(12):
        nodeRho = np.bincount(cell.flat, weights=np.repeat(cellRho*area, 3), minlength=NN)
        nodeRho /= sumArea
        cellRho = np.sum(nodeRho[cell], axis=-1)/3 
    cellRho = (cellRho/np.max(cellRho))**(gamma)
    return cellRho, normal