Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

外齿轮齿面内法向计算 #1449

Merged
merged 5 commits into from
Jan 14, 2025
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
gearx: 修复bug,齿面节点获取函数
  • Loading branch information
concha-k-chen committed Jan 14, 2025
commit aa6d57b88848bb40605524dff455bb5d866d70fd
144 changes: 81 additions & 63 deletions app/gearx/gear.py
Original file line number Diff line number Diff line change
@@ -62,7 +62,7 @@ def __init__(self, m_n, z, alpha_n, beta, x_n, hac, cc, rcc, jn, n1, n2, n3, na,
self.mesh = None
self.hex_mesh = None
self.target_quad_mesh = None
self.targer_hex_mesh = None
self.target_hex_mesh = None
self.gear_type = 0
self.is_max_angle = False
self.single_node_num = 0
@@ -323,45 +323,70 @@ def get_profile_node_index(self, tooth_tag=None):
tooth_tag = list(tooth_tag)
else:
raise TypeError('The tooth_tag must be an integer or a list of integers.')
quad_mesh = self.mesh
hex_mesh = self.hex_mesh
# 构建齿面参考法向
rotation = self.rotation_direction*self.beta
tooth_tangnet = np.array([0, np.sin(rotation), np.cos(rotation)])
tooth_normal_right = np.array([0, np.cos(rotation), -np.sin(rotation)])
tooth_normal_left = np.array([0, -np.cos(rotation), np.sin(rotation)])

quad_mesh = self.target_quad_mesh
hex_mesh = self.target_hex_mesh
t_NN = quad_mesh.number_of_nodes()
t_NC = quad_mesh.number_of_cells()
n1 = self.n1
nw = self.nw
hex_node = hex_mesh.node

is_bd_cell = quad_mesh.boundary_cell_flag()
domain_flag_left = quad_mesh.celldata["cell_domain_tag"] == 6
domain_flag_right = quad_mesh.celldata["cell_domain_tag"] == 5
domain_flag_right = quad_mesh.celldata["cell_domain_tag"] == 6
domain_flag_left = quad_mesh.celldata["cell_domain_tag"] == 5
if tooth_tag is None:
folder = self.z
cell_flag_left = is_bd_cell & domain_flag_left
cell_flag_right = is_bd_cell & domain_flag_right
cell_flag_left = is_bd_cell & domain_flag_left
else:
folder = len(tooth_tag)
tooth_flag = np.isin(quad_mesh.celldata["cell_tooth_tag"], tooth_tag)

cell_flag_left = is_bd_cell & domain_flag_left & tooth_flag
cell_flag_right = is_bd_cell & domain_flag_right & tooth_flag

# 右侧齿面
# TODO: 考虑是否给单元局部起始节点编号加检测,即不一定从 0 和 2 开始
cell_idx_right = np.where(cell_flag_right)[0].reshape(folder, -1)[..., 0:n1]
tooth_profile_cell_right = quad_mesh.cell[cell_idx_right]
tooth_profile_node_right = np.zeros((folder, nw + 1, n1 + 1), dtype=np.int32)
tooth_profile_node_right[..., 0, 0:n1] = tooth_profile_cell_right[..., :, 0]
tooth_profile_node_right[..., 0, -1] = tooth_profile_cell_right[..., -1, 1]
cell_flag_left = is_bd_cell & domain_flag_left & tooth_flag

# 左侧齿面
cell_idx_left = np.flip(np.flip(np.where(cell_flag_left)[0]).reshape(folder, -1)[..., 0:n1], axis=0)
# TODO: 考虑是否给单元局部起始节点编号加检测,即不一定从 0 和 2 开始
cell_idx_left = np.where(cell_flag_left)[0].reshape(folder, -1)[..., 0:n1]
tooth_profile_cell_left = quad_mesh.cell[cell_idx_left]
tooth_profile_node_left = np.zeros((folder, nw + 1, n1 + 1), dtype=np.int32)
tooth_profile_node_left[..., 0, 0:n1] = tooth_profile_cell_left[..., :, 2]
tooth_profile_node_left[..., 0, 0:n1] = tooth_profile_cell_left[..., :, 0]
tooth_profile_node_left[..., 0, -1] = tooth_profile_cell_left[..., -1, 1]
# 左侧齿面节点所在单元记录
left_cell_idx = np.zeros_like(tooth_profile_node_left)
left_cell_idx[..., 0, 0:n1] = cell_idx_left[..., :]
left_cell_idx[..., 0, -1] = cell_idx_left[..., -1]

# 右侧齿面
cell_idx_right = np.flip(np.flip(np.where(cell_flag_right)[0]).reshape(folder, -1)[..., 0:n1], axis=0)
tooth_profile_cell_right = quad_mesh.cell[cell_idx_right]
tooth_profile_node_right = np.zeros((folder, nw + 1, n1 + 1), dtype=np.int32)
tooth_profile_node_right[..., 0, 0:n1] = tooth_profile_cell_right[..., :, 2]
tooth_profile_node_right[..., 0, -1] = tooth_profile_cell_right[..., -1, 1]
# 右侧齿面节点所在单元记录
right_cell_idx = np.zeros_like(tooth_profile_node_right)
right_cell_idx[..., 0, 0:n1] = cell_idx_right[..., :]
right_cell_idx[..., 0, -1] = cell_idx_right[..., -1]

for i in range(1, nw + 1):
tooth_profile_node_right[..., i, :] = tooth_profile_node_right[..., 0, :] + i * t_NN
tooth_profile_node_left[..., i, :] = tooth_profile_node_left[..., 0, :] + i * t_NN
tooth_profile_node_right[..., i, :] = tooth_profile_node_right[..., 0, :] + i * t_NN
for i in range(1, nw):
left_cell_idx[..., i, :] = left_cell_idx[..., 0, :] + i * t_NC
right_cell_idx[..., i, :] = right_cell_idx[..., 0, :] + i * t_NC
left_cell_idx[..., -1, :] = left_cell_idx[..., 0, :] + (nw-1) * t_NC
right_cell_idx[..., -1, :] = right_cell_idx[..., 0, :] + (nw-1) * t_NC

# 计算齿廓节点内法向
cell = hex_mesh.cell
right_profile_node_face_normal = get_face_normal_with_reference(hex_mesh[right_cell_idx], tooth_normal_right)
left_profile_node_face_normal = get_face_normal_with_reference(hex_mesh[left_cell_idx], tooth_normal_left)

return (tooth_profile_node_right, tooth_profile_node_left), (
hex_node[tooth_profile_node_right], hex_node[tooth_profile_node_left])
@@ -388,39 +413,62 @@ def set_target_tooth(self, target_tooth_tag, get_wheel=False):
total_quad_mesh = self.mesh
total_hex_mesh = self.hex_mesh

if get_wheel:
total_node = total_quad_mesh.node
total_cell = total_quad_mesh.cell
total_number_of_nodes = total_quad_mesh.number_of_nodes()
total_number_of_cells = total_quad_mesh.number_of_cells()
cell_tooth_tag = total_quad_mesh.celldata["cell_tooth_tag"]
else:
total_node = total_hex_mesh.node
total_cell = total_hex_mesh.cell
total_number_of_nodes = total_hex_mesh.number_of_nodes()
total_number_of_cells = total_hex_mesh.number_of_cells()
cell_tooth_tag = total_hex_mesh.celldata["cell_tooth_tag"]
# 构建目标齿端面网格
total_node = total_quad_mesh.node
total_cell = total_quad_mesh.cell
total_number_of_nodes = total_quad_mesh.number_of_nodes()
total_number_of_cells = total_quad_mesh.number_of_cells()
cell_tooth_tag = total_quad_mesh.celldata["cell_tooth_tag"]
cell_domain_tag = total_quad_mesh.celldata["cell_domain_tag"]
# 获取目标齿相关单元
is_target_cell = np.zeros(total_number_of_cells, dtype=bool)
for tag in target_tooth_tag:
is_target_cell |= (cell_tooth_tag == tag)
target_cell_idx = np.where(is_target_cell)[0]
target_cell = total_cell[target_cell_idx]

# 标记目标齿面相关节点
is_target_node = np.zeros(total_number_of_nodes, dtype=bool)
is_target_node[target_cell] = True
target_node_idx = np.where(is_target_node)[0]
target_node = total_node[target_node_idx]

# 构建节点映射
node_idx_map = np.zeros(total_number_of_nodes, dtype=int)
node_idx_map[target_node_idx] = np.arange(len(target_node_idx))
# 单元映射
target_cell = node_idx_map[target_cell]
self.target_quad_mesh = QuadrangleMesh(target_node, target_cell)
target_cell_domain_tag = cell_domain_tag.reshape(self.z, -1)[target_tooth_tag]
target_cell_tooth_tag = cell_tooth_tag.reshape(self.z, -1)[target_tooth_tag]
self.target_quad_mesh.celldata["cell_domain_tag"] = target_cell_domain_tag.reshape(-1)
self.target_quad_mesh.celldata["cell_tooth_tag"] = target_cell_tooth_tag.reshape(-1)

if not get_wheel:
self.target_hex_mesh = HexahedronMesh(target_node, target_cell)
node = target_node
cell = target_cell
beta = self.beta
r = self.r
tooth_width = self.tooth_width
nw = self.nw
rotation_direction = self.rotation_direction
# 数据处理,将二维点转换为三维点
new_node = np.zeros((len(node), 3))
new_node[:, 0:2] = node
one_section_node_num = len(new_node)

# 创建齿轮整体网格
# 拉伸节点
volume_node = sweep_points(new_node, beta, r, tooth_width, nw, rotation_direction).reshape(-1, 3)
# 将端面四边形单元拉伸为六面体单元
volume_cell = np.zeros((nw, len(cell), 8), dtype=np.int64)
cell_domain_tag = np.zeros((nw, len(cell)))
cell_tooth_tag = np.zeros((nw, len(cell)))
# 填充单元的节点索引
for i in range(nw):
volume_cell[i, :, 0:4] = cell + i * one_section_node_num
volume_cell[i, :, 4:8] = cell + (i + 1) * one_section_node_num
volume_cell = volume_cell.reshape(-1, 8)

self.target_hex_mesh = HexahedronMesh(volume_node, volume_cell)
else:
df = self.d_f
total_tooth_node_num = len(target_node)
@@ -824,11 +872,6 @@ def generate_mesh(self):
z = self.z
# 获取齿廓与过渡曲线点列
points = self.get_profile_points()
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 10))
plt.plot(points[:51, 0], points[:51, 1], 'o-', label='Profile Points')
plt.axis('equal')
plt.show()
r_inner = self.inner_diam / 2

one_tooth_angle = abs(delta_angle_calculator(points[0, :2], points[n1 + n2 + 1, :2], input_type="vector"))
@@ -961,16 +1004,6 @@ def generate_mesh(self):
np.linspace(key_points[edge[25, 0]], key_points[edge[25, 1]], na2 + 1),
np.linspace(key_points[edge[26, 0]], key_points[edge[26, 1]], na1 + 1)
]
import matplotlib.pyplot as plt
# 验证 line
fig, ax = plt.subplots()
for l in line:
ax.plot(l[:, 0], l[:, 1])
# l = line[9]
# ax.plot(l[:, 0], l[:, 1])
plt.axis('equal')
plt.axis('off')
plt.show()

# 构建子区域半边数据结构
# boundary_edge = np.array([0, 1, 4, 5, 6, 7, 8, 9, 10, 11, 3, 2, 12, 13])
@@ -2071,11 +2104,6 @@ def generate_mesh(self):
t_mesh.celldata['cell_tooth_tag'] = cell_tooth_tag
else:
points = self.get_profile_points()
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 10))
plt.plot(points[:31, 0], points[:31, 1], 'o-', label='Profile Points')
plt.axis('equal')
plt.show()
# 构建关键点
angle_kp0 = pi / 2 - 2 * pi / z / 2
angle_kp4 = pi / 2 + 2 * pi / z / 2
@@ -2233,16 +2261,6 @@ def generate_mesh(self):
r_kp2 * np.concatenate([cos(delta_kp23_17), sin(delta_kp23_17)], axis=1),
r_kp2 * np.concatenate([cos(delta_kp17_20), sin(delta_kp17_20)], axis=1)
]
import matplotlib.pyplot as plt
# 验证 line
fig, ax = plt.subplots()
for l in line:
ax.plot(l[:, 0], l[:, 1])
# l = line[9]
# ax.plot(l[:, 0], l[:, 1])
plt.axis('equal')
plt.axis('off')
plt.show()

# boundary_edge = np.array([0, 1, 2, 6, 30, 31, 7, 5, 4, 3, 8, 9, 10, 11, 28, 29, 12, 13, 14, 15])
# half_edge = subdomain_divider(line, key_points, edge, boundary_edge)
60 changes: 56 additions & 4 deletions app/gearx/test/test_gear_system.py
Original file line number Diff line number Diff line change
@@ -35,6 +35,14 @@ def test_external_gear(self):
nw = data['nw']
tooth_width = data['tooth_width']
inner_diam = data['inner_diam'] # 轮缘内径
# m_n = 2.25
# z = 64
# alpha_n = 17.5
# beta = 30
# x_n = 0.4
# hac = 1
# cc = 1.4
# inner_diam = 10
chamfer_dia = data['chamfer_dia'] # 倒角高度(直径)

external_gear = ExternalGear(m_n, z, alpha_n, beta, x_n, hac, cc, rcc, jn, n1, n2, n3, na, nf, nw, chamfer_dia,
@@ -43,7 +51,7 @@ def test_external_gear(self):
quad_mesh = external_gear.generate_mesh()
# quad_mesh.to_vtk(fname='../data/external_quad_mesh.vtu')
hex_mesh = external_gear.generate_hexahedron_mesh()
hex_mesh.to_vtk(fname='../data/test_param_external_hex_mesh.vtu')
# hex_mesh.to_vtk(fname='../data/test_param_external_hex_mesh.vtu')
#
# with open('external_gear.pkl', 'wb') as f:
# pickle.dump({'quad_mesh': quad_mesh, 'gear': external_gear, 'hex_mesh': hex_mesh}, f)
@@ -54,8 +62,8 @@ def test_external_gear(self):
# quad_mesh_from_cpp = QuadrangleMesh(node_from_cpp, cell_from_cpp)
# quad_mesh_from_cpp.to_vtk(fname='external_quad_mesh_cpp.vtu')

# with open('../data/external_gear.pkl', 'wb') as f:
# pickle.dump({'external_gear': external_gear, 'hex_mesh': hex_mesh, 'quad_mesh': quad_mesh}, f)
with open('../data/external_gear_new.pkl', 'wb') as f:
pickle.dump({'external_gear': external_gear, 'hex_mesh': hex_mesh, 'quad_mesh': quad_mesh}, f)

def test_internal_gear(self):
with open('../data/internal_gear_data.json', 'r') as file:
@@ -83,6 +91,9 @@ def test_internal_gear(self):
internal_gear = InternalGear(m_n, z, alpha_n, beta, x_n, hac, cc, rcc, jn, n1, n2, n3, na, nf, nw, outer_diam, z_cutter,
xn_cutter, tooth_width)

hex_mesh = internal_gear.generate_hexahedron_mesh()
hex_mesh.to_vtk(fname='../data/test_param_internal_hex_mesh.vtu')

t = np.array([0.1, 0.2, 25])
p1 = internal_gear.get_involute_points(t)
p1_dis = internal_gear.get_tip_intersection_points(t)
@@ -98,7 +109,7 @@ def test_internal_gear(self):

# p = internal_gear.get_profile_points()

quad_mesh = internal_gear.generate_mesh()
# quad_mesh = internal_gear.generate_mesh()
# quad_mesh.to_vtk(fname='../data/internal_quad_mesh.vtu')
# r = internal_gear.r
# # hex_mesh = generate_hexahedral_mesh(quad_mesh, internal_gear.beta, r, tooth_width, nw)
@@ -480,6 +491,47 @@ def test_get_part_external_gear(self):
part_external_gear_mesh.to_vtk(fname='../data/part_external_gear_mesh.vtu')
part_external_gear_mesh_wheel.to_vtk(fname='../data/part_external_gear_hex_mesh_wheel.vtu')

def test_get_profile_face_normal(self):
with open('../data/external_gear_new.pkl', 'rb') as f:
data = pickle.load(f)
external_gear = data['external_gear']
hex_mesh = data['hex_mesh']
quad_mesh = data['quad_mesh']

mesh = external_gear.set_target_tooth([0, 1, 18])

GD = mesh.geo_dimension()
NC = mesh.number_of_cells()
print(f"NC: {NC}")
NN = mesh.number_of_nodes()
print(f"NN: {NN}")
node = mesh.entity('node')
cell = mesh.entity('cell')
print(f"cell: {cell}")
print(f"node: {node}")

# 齿面上的节点索引和坐标
node_indices_tuple, noe_coord_tuple = external_gear.get_profile_node_index(tooth_tag=0)
node_indices_right = node_indices_tuple[0].reshape(-1, 1)
node_indices_left = node_indices_tuple[1].reshape(-1, 1)
node_indices = np.concatenate([node_indices_right, node_indices_left], axis=0) # (NPN, 1)
node_coord_right = noe_coord_tuple[0].reshape(-1, 3)
node_coord_left = noe_coord_tuple[1].reshape(-1, 3)
node_coord = np.concatenate([node_coord_right, node_coord_left], axis=0) # (NPN, GD)
# 齿面上的节点数
NPN = node_indices.shape[0]
# TODO 齿面上节点的内法线方向
face_normal = np.zeros((NPN, 3), np.float64)
for i, t_node in enumerate(node_coord):
t1, face_normal[i], t2 = external_gear.find_node_location_kd_tree(t_node)
average_normal = np.mean(face_normal, axis=0)
average_normal /= np.linalg.norm(average_normal)
threshold = 0.1
for i in range(len(face_normal)):
deviation = np.linalg.norm(face_normal[i] - average_normal)
if deviation > threshold:
face_normal[i] = average_normal




28 changes: 28 additions & 0 deletions app/gearx/utils.py
Original file line number Diff line number Diff line change
@@ -717,3 +717,31 @@ def face_normal_bilinear(cell_nodes, u, v, w, error=1e-3, is_normalize=True):
return normal
else:
return np.cross(Pu, Pv)

# 给定参考法向计算某点的法向
def get_face_normal_with_reference(cell, reference_normal, error=1e-3, is_normalize=True):
'''
给定一组六面体单元,以及一个参考法向,计算该单元的六个面的法向,并返回与参考法向最接近的法向
:param cell: 计算单元
:param reference_normal: 参考法向
:param error: 误差限制
:param is_normalize: 是否归一化
'''
NC = cell.shape[0]
face_normal = np.zeros((NC, 3))
parameters = np.array([
[0.5, 0.5, 0],
[0.5, 0, 0.5],
[0, 0.5, 0.5],
[0.5, 0.5, 1],
[0.5, 1, 0.5],
[1, 0.5, 0.5]
])
for i in range(NC):
for j in range(6):
temp_normal = face_normal_bilinear(cell[i], parameters[j, 0], parameters[j, 1], parameters[j, 2], error, is_normalize)
vec_err = np.sqrt(np.sum((temp_normal.reshape(-1) - reference_normal)**2))
if vec_err < error:
face_normal[i] = temp_normal
break
return face_normal