计算由 scipy 的 Delaunay 函数生成的四面体网格的雅可比矩阵

如何解决计算由 scipy 的 Delaunay 函数生成的四面体网格的雅可比矩阵

我正在尝试使用 Delaynay 的函数 scipy 生成四面体网格。从提供的源代码 here 中，我做了以下内容：

import math
import random
import numpy as np
from scipy.spatial import delaunay

# geometry
height = 2.0
depth  = 2.0
width  = 2.0

# random transformations of the original mesh
random.seed(1234)
variation = .00 # as a decimal for a percent

h = 0.4

x_layer = int(width / h)+1
y_layer = int(depth / h)+1
z_layer = int(height / h)+1

points = np.zeros([x_layer*y_layer*z_layer,3])

count = 0
for orig_x in range(x_layer):
    for orig_y in range(y_layer):
        for orig_z in range(z_layer):

            rand_pert_x = random.uniform(-variation,variation)
            rand_pert_y = random.uniform(-variation,variation)
            rand_pert_z = random.uniform(-variation,variation)
    
            x = orig_x*h
            y = orig_y*h
            z = orig_z*h
     
            if (x==0 or orig_x==x_layer-1) or (y==0 or orig_y==y_layer-1) or (z==0 or orig_z==z_layer-1):
                points[count,:] = np.array([x,y,z])
            else:
                points[count,:] = np.array([x+h*rand_pert_x,y+h*rand_pert_y,z+h*rand_pert_z])
            count += 1

tet = delaunay(points,qhull_options="Qt")

即计划是让点随机扰动以创建低质量的网格。目前，variation 设置为 0.0，因此代码将生成 3D 矩形网格。

我想测试网格的质量，所以我要看的一件事是每个元素的雅可比矩阵。有没有办法从 delaunay 函数的输出中计算出来？

解决方法

以下代码计算数组

假设一个四面体的点为 pj=(xj,yj,zj) 为 j=0,1,2,3，与之对应的雅可比矩阵为（参见示例 here）：

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

以下函数返回与每个四面体和每个矩阵的行列式值对应的雅可比矩阵数组。

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,#  [y1-y0,#  [z1-z0,z3-z0]]
    # 
    a = simp_pts[:,0] - simp_pts[:,0]  # x1-x0
    b = simp_pts[:,1] - simp_pts[:,1]  # y1-y0
    c = simp_pts[:,2] - simp_pts[:,2]  # z1-z0
    d = simp_pts[:,0]  # x2-x0
    e = simp_pts[:,1]  # y2-y0
    f = simp_pts[:,2]  # z2-z0
    g = simp_pts[:,0]  # x3-x0
    h = simp_pts[:,1]  # y3-y0
    i = simp_pts[:,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))
    jacobian_array[:,0] = a
    jacobian_array[:,1] = b
    jacobian_array[:,2] = c
    jacobian_array[:,0] = d
    jacobian_array[:,1] = e
    jacobian_array[:,2] = f
    jacobian_array[:,0] = g
    jacobian_array[:,1] = h
    jacobian_array[:,2] = i

    return jacobian_array,determinants