How can I avoid tetrahedra lying in a plane using scipy.Delaunay?

I'm trying to create a tetrahedral mesh in Python 3.7.3.. But some tetrahedra are flat, i.e. their vertices are on a flat surface.

import numpy as np
from scipy.spatial import Delaunay

# Coordinates of 3x3x3 equally distant points of a cube
x = np.linspace(0, 1, 3)
X1, X2, X3 = np.meshgrid(x,x,x)
vertices = np.hstack([X1.reshape(-1,1),X2.reshape(-1,1),X3.reshape(-1,1)])

# Using Delaunay
tri = Delaunay(vertices).simplices

# tetrahedra
simplices = vertices[tri,:]

As you can see below, already with the third tetrahedron all y-coordinates are 0.5. This leads later to a singular matrix.

print(simplices[0:3])
[[[1.  0.5 1. ]
  [0.5 1.  0.5]
  [1.  0.5 0.5]
  [0.5 0.5 0.5]]

 [[1.  0.5 1. ]
  [1.  1.  0.5]
  [0.5 1.  0.5]
  [1.  0.5 0.5]]

 [[1.  0.5 1. ]
  [0.5 0.5 1. ]
  [1.  0.5 0.5]
  [0.5 0.5 0.5]]]

Do you know how I can work around this problem? Thank you very much.

Solution

Unlike 2D, for 3D there is no known algorithm that definitely generates a Delaunay triangulation of a given domain. There are however some mesh generation packages which produce pretty good tetrahedral meshes. For example

(Disclaimer: I'm the author of pygmsh and pygalmesh.)