pyvistaによる3Dデータ同士の足し引きについて

前提

pythonのpyvistaを用いて、２つの３Dデータ同士の足し引きを行おうと考えています。

実現したいこと

以下のような平行六面体の中から、

固定位置にある球体をくりぬいていき、

最終的に平行六面体の中にクレーターができたような立体を作りたいと考えています。

発生している問題・エラーメッセージ

しかし、以下のようなエラーが表示され、思うように立体が出力されません。

ERROR:root:No points/cells to operate on
Traceback (most recent call last):
File "C:\Users*\cad3.py", line 34, in <module>
subtracted=subtracted-pv.Sphere(radius=1.5, center=(x[0], x[1], x[2]))
File "C:\Users*\pyvista\core\pointset.py", line 928, in sub
return self.boolean_difference(cutting_mesh)
File "C:\Users*\pyvista\core\filters\poly_data.py", line 316, in boolean_difference
return self._boolean('difference', other_mesh, tolerance, progress_bar=progress_bar)
File "C:\Users*\pyvista\core\filters\poly_data.py", line 78, in _boolean
raise NotAllTrianglesError("Make sure both the input and output are triangulated.")
pyvista.core.errors.NotAllTrianglesError: Make sure both the input and output are triangulated.

該当のソースコード

ソースコードは以下の通りとなります。

import numpy as np
import pyvista as pv
from pyvista import examples

#平行六面体
a=[10.3559103, 0, 0]
b=[5.896951249, 8.512980919, 0]
c=[5.896951249, 3.088725841, 7.932881999]
vertices=np.array([[0, 0, 0], a, [x+y for x, y in zip(a, b)], b, c, [x+y for x, y in zip(a, c)], [x+y+z for x, y, z in zip(a, b, c)], [x+y for x, y in zip(b, c)]])
faces=np.hstack([[4, 0, 1, 2, 3], [4, 4, 5, 6, 7], [4, 0, 1, 5, 4], [4, 3, 2, 6, 7], [4, 0, 4, 7, 3], [4, 1, 5, 6, 2]])
hexahedron=pv.PolyData(vertices, faces).triangulate()

#球体の座標
particle_coordinates=[
[6.09692140007496, 7.11117400736113, 0.90998651328237],
[8.81611700568463, 4.18411507652898, 4.41893502931168],
[12.5421432054935, 4.99850518084874, 7.46484920848187],
[5.7222559705099, 1.92865547036078, 2.66685804719921],
[15.5954635493239, 4.07668654519134, 7.87222596090579],
[13.0365970393955, 7.2950214675815, 4.89336477824595],
[12.1461765720701, 7.69670673623446, 7.6901335880258],
[14.9977694749948, 8.27629391793453, 7.31868302287639],
[13.6403367137098, 8.26192585787121, 1.63089634553616],
[14.4373121118686, 4.67871533316406, 5.11464550601519],
]

#平行六面体から球体を引いていく処理
subtracted=hexahedron
for x in particle_coordinates:
subtracted=subtracted-pv.Sphere(radius=1.5, center=(x[0], x[1], x[2]))

#立体の表示
pl=pv.Plotter()
pl.add_mesh(subtracted)
pl.show()

おそらく初歩的な内容かと存じますが、どなたか原因のわかる方がいらしたら、お教えいただけますでしょうか。