Use skfem to solve Laplace eigenvalue problem on triangular mesh on 3D surface

Question

I am trying to build a toolchain to solve the Laplace eigenvalue problem on a three-dimensional surface. I first define a point cloud. Then I generate a surface from the grid points using pyvista. Then I bring the points and faces into the correct form before trying to set up the finite element approximation to solve the PDE.

# Import packages
import numpy as np
import pyvista as pv
from math import pi
import skfem as fem

# Transform faces into right structure
def clean_faces(faces):
    # Split arrays
    subarrays = np.array_split(faces, len(faces) // 4)

    # Transform to array
    subarrays_list = [subarray.tolist() for subarray in subarrays]

    # Pick only the points that define the triangles
    arr = [sublist[1:] for sublist in subarrays_list]

    # Transform back to array
    cells = np.array(arr)
    return cells

# Clean up data
def _sanitize(points, cells):
    uvertices, uidx = np.unique(cells, return_inverse=True)
    cells = uidx.reshape(cells.shape)
    points = points[uvertices]
    return points, cells

# Intialize parameters
r = 1;
nn = 20;
GridPoints = []


# Define point cloud of a sphere
for jj in range(nn):
    for ii in range(nn):
        px = r * np.sin(pi*ii/nn) * np.cos(2*pi*jj/nn)
        py = r * np.sin(pi*ii/nn) * np.sin(2*pi*jj/nn)
        pz = r * np.cos(pi*ii/nn)
        GridPoints.append([px,py,pz])

# Transform GridPoints into array
points = np.array(GridPoints)

# Set up array as pyvist point cloud
cloud = pv.PolyData(points)

# Reconstruct surface of the gridpoints
surf = cloud.reconstruct_surface()

# Extract and clean surfaces triangles
tri = clean_faces(surf.faces)

# Sanitize data
points, cells = _sanitize(surf.points, tri)

# Bring inputs in the right form
p = np.array([points[:,0],  points[:,1],  points[:,2]], dtype=np.float64)
t = np.array(cells, dtype=np.float64).T

# Create fem mesh
myMesh = fem.MeshTri(p,t).refined(2)

# piecewise linear finite elements
e = fem.ElementTriP1()

# create mapping for the finite element approximation
basis = fem.CellBasis(myMesh, e)

When I want to create a mapping for the finite element approximation I get an error. I get the following error message.

> IndexError: index 3 is out of bounds for axis 0 with size 3

I tried using different fem elements like tetrahedral elements (fem.ElementTetP1(), fem.MeshTet(p, t), etc.) with no success. I also tried to clean up the data as much as possible. For example by excluding the points that are not connected to rest to the triangular mesh. Also with no success.