计算由 scipy 的 Delaunay 函数生成的四面体网格的雅可比矩阵
Calculate the Jacobian of the tetrahedral mesh generated from scipy's Delaunay function
我正在尝试使用 scipy
的函数 Delaynay
生成四面体网格。根据提供的源代码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]]
下面的函数returns每个四面体对应的雅可比矩阵数组和每个矩阵的行列式值
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
我正在尝试使用 scipy
的函数 Delaynay
生成四面体网格。根据提供的源代码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]]
下面的函数returns每个四面体对应的雅可比矩阵数组和每个矩阵的行列式值
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