Calculate The Jacobian Of The Tetrahedral Mesh Generated From Scipy's Delaunay Function

I am trying to use the function Delaynay of scipy to generate a tetrahedral mesh. From the source code provided here, I have make something as followed: import math import random i

Solution 1:

The following code computes an array

Assuming a tetrahedron with points pj=(xj, yj, zj) for j=0, 1, 2, 3, the Jacobian matrix corresponding to it is (see for example here):

   [[x1-x0, x2-x0, x3-x0], 
J = [y1-y0, y2-y0, y3-y0], 
    [z1-z0, z2-z0, z3-z0]]

The following function returns an array of Jacobian matrices corresponding to each of the tetrahedra and the determinant value of each matrix.

def compute_delaunay_tetra_jacobians(dt):
    """
    Compute the Jacobian matrix and its determinant for each tetrahedron in the Delaunay triangulation.
    :param dt: the Delaunay triangulation
    :return: array of shape (n, 3, 3) of jacobian matrices such that jacoboian_array[i, :, :] is the 3x3 Jacobian matrix
             array of n values of the determinant of the jacobian matrix
    """
    simp_pts = dt.points[dt.simplices]
    # (n, 4, 3) array of tetrahedra points where simp_pts[i, j, k] holds the k'th coordinate (x, y or z)
    # of the j'th 3D point (of four) of the i'th tetrahedron
    assert simp_pts.shape[1] == 4 and simp_pts.shape[2] == 3

    # building the 3x3 jacobian matrix with entries:
    # [[x1-x0, x2-x0, x3-x0], 
    #  [y1-y0, y2-y0, y3-y0], 
    #  [z1-z0, z2-z0, z3-z0]]
    # 
    a = simp_pts[:, 1, 0] - simp_pts[:, 0, 0]  # x1-x0
    b = simp_pts[:, 1, 1] - simp_pts[:, 0, 1]  # y1-y0
    c = simp_pts[:, 1, 2] - simp_pts[:, 0, 2]  # z1-z0
    d = simp_pts[:, 2, 0] - simp_pts[:, 0, 0]  # x2-x0
    e = simp_pts[:, 2, 1] - simp_pts[:, 0, 1]  # y2-y0
    f = simp_pts[:, 2, 2] - simp_pts[:, 0, 2]  # z2-z0
    g = simp_pts[:, 3, 0] - simp_pts[:, 0, 0]  # x3-x0
    h = simp_pts[:, 3, 1] - simp_pts[:, 0, 1]  # y3-y0
    i = simp_pts[:, 3, 2] - simp_pts[:, 0, 2]  # z3-z0

    determinants = a*(e*i - f*h) - b*(d*i - f*g) + c*(d*h - e*g)

    n = simp_pts.shape[0]
    jacobian_array = np.empty((n, 3, 3))
    jacobian_array[:, 0, 0] = a
    jacobian_array[:, 0, 1] = b
    jacobian_array[:, 0, 2] = c
    jacobian_array[:, 1, 0] = d
    jacobian_array[:, 1, 1] = e
    jacobian_array[:, 1, 2] = f
    jacobian_array[:, 2, 0] = g
    jacobian_array[:, 2, 1] = h
    jacobian_array[:, 2, 2] = i

    return jacobian_array, determinants