def METHOD_NAME(meshes, method: str = "cot"):
    """
    Naive implementation of laplacian smoothing wit cotangent weights.
    """
    verts_packed = meshes.verts_packed()  # (sum(V_n), 3)
    faces_packed = meshes.faces_packed()  # (sum(F_n), 3)
    V = verts_packed.shape[0]
    L = torch.zeros((V, V), dtype=torch.float32, device=meshes.device)
    inv_areas = torch.zeros((V, 1), dtype=torch.float32, device=meshes.device)
    for f in faces_packed:
        v0 = verts_packed[f[0], :]
        v1 = verts_packed[f[1], :]
        v2 = verts_packed[f[2], :]
        A = (v1 - v2).norm()
        B = (v0 - v2).norm()
        C = (v0 - v1).norm()
        s = 0.5 * (A + B + C)
        face_area = (s * (s - A) * (s - B) * (s - C)).clamp_(min=1e-12).sqrt()
        inv_areas[f[0]] += face_area
        inv_areas[f[1]] += face_area
        inv_areas[f[2]] += face_area
        A2, B2, C2 = A * A, B * B, C * C
        cota = (B2 + C2 - A2) / face_area / 4.0
        cotb = (A2 + C2 - B2) / face_area / 4.0
        cotc = (A2 + B2 - C2) / face_area / 4.0
        L[f[1], f[2]] += cota
        L[f[2], f[0]] += cotb
        L[f[0], f[1]] += cotc
        # symetric
        L[f[2], f[1]] += cota
        L[f[0], f[2]] += cotb
        L[f[1], f[0]] += cotc
    idx = inv_areas > 0
    inv_areas[idx] = 1.0 / inv_areas[idx]
    norm_w = L.sum(dim=1, keepdims=True)
    L_sum = norm_w.clone()
    idx = norm_w > 0
    norm_w[idx] = 1.0 / norm_w[idx]
    if method == "cotcurv":
        loss = (L.mm(verts_packed) - L_sum * verts_packed) * inv_areas * 0.25
        loss = loss.norm(dim=1)
    else:
        loss = L.mm(verts_packed) * norm_w - verts_packed
        loss = loss.norm(dim=1)
    weights = torch.zeros(V, dtype=torch.float32, device=meshes.device)
    for v in range(V):
        weights[v] = meshes.num_verts_per_mesh()[
            meshes.verts_packed_to_mesh_idx()[v]
        ]
    weights = 1.0 / weights
    loss = loss * weights
    return loss.sum() / len(meshes)