Merge branch 'develop' of github.com:weihuayi/fealpy into develop

weihuayi · Jan 17, 2025 · c5c461a · c5c461a
2 parents 8140707 + 7d4e36f
commit c5c461a
Show file tree

Hide file tree

Showing 17 changed files with 651 additions and 191 deletions.
diff --git a/app/gearx/gear.py b/app/gearx/gear.py
@@ -595,25 +595,26 @@ def set_target_tooth(self, target_tooth_tag, get_wheel=False):
                 if (target_tooth_tag[0] == 0 and (target_tooth_tag[-1] != (z - 1))) or target_tooth_tag[0] != 0:
                     node_start = target_node[-single_node_num]
                     node_end = target_node[0]
-                    node_temp1 = target_node[1]
+                    # node_temp1 = target_node[1]
                     node_temp2 = target_node[2]
                 else:
                     node_start = target_node[start_num * single_node_num + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)]
                     node_end = target_node[end_num * single_node_num + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)]
-                    node_temp1 = target_node[end_num * single_node_num + 1  + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)]
+                    # node_temp1 = target_node[end_num * single_node_num + 1 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)]
                     node_temp2 = target_node[end_num * single_node_num + 2 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)]
                 angle_start = np.arctan2(node_start[1], node_start[0]) % (2 * pi)
                 angle_end = np.arctan2(node_end[1], node_end[0]) % (2 * pi)
-                r_temp1 = np.sqrt(node_temp1[0] ** 2 + node_temp1[1] ** 2)
+                # r_temp1 = np.sqrt(node_temp1[0] ** 2 + node_temp1[1] ** 2)
                 r_temp2 = np.sqrt(node_temp2[0] ** 2 + node_temp2[1] ** 2)
 
                 if angle_end < angle_start:
                     angle_end += 2 * pi
                 phi = np.linspace(angle_start, angle_end, 2 * na2 * (self.z - tooth_num))[1:-1]
-                r1 = np.linspace(self.r_a, r_temp1, n1 + 1)
-                r2 = np.linspace(r_temp1, r_temp2, n2 + 1)
-                r3 = np.linspace(r_temp2, self.outer_diam/2, n3 + 1)
-                r = np.concatenate((r1, r2[1:], r3[1:]))
+                # r1 = np.linspace(self.r_a, r_temp1, n1 + 1)
+                # r2 = np.linspace(r_temp1, r_temp2, n2 + 1)
+                # r3 = np.linspace(r_temp2, self.outer_diam/2, n3 + 1)
+                r = np.linspace(r_temp2, self.outer_diam/2, n3 + 1)
+                # r = np.concatenate((r1, r2[1:], r3[1:]))
 
                 x = np.einsum('r,p->pr', r, np.cos(phi))[..., None]
                 y = np.einsum('r,p->pr', r, np.sin(phi))[..., None]
@@ -623,52 +624,52 @@ def set_target_tooth(self, target_tooth_tag, get_wheel=False):
                 wheel_node_idx = np.zeros((len(phi) + 2, len(r)), dtype=np.int64)
                 wheel_node_idx[1:-1, :] = np.arange(total_tooth_node_num,
                                                     total_tooth_node_num + len(wheel_node)).reshape(-1, len(r))
-                wheel_node_idx[0, 0] = start_num * single_node_num + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
-                wheel_node_idx[0, 1:n1] = np.arange(
-                    start_num * single_node_num + number_of_key_points + (n1 - 1) + (n2 - 1) + (n3 - 1),
-                    start_num * single_node_num + number_of_key_points + 2 * (n1 - 1) + (n2 - 1) + (n3 - 1))[::-1]
-                wheel_node_idx[0, n1] = start_num * single_node_num + 1 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
-                wheel_node_idx[0, n1 + 1:n1 + n2] = np.arange(
-                    start_num * single_node_num + number_of_key_points + 2 * (n1 - 1) + (n2 - 1) + (n3 - 1),
-                    start_num * single_node_num + number_of_key_points + 2 * (n1 - 1) + 2 * (n2 - 1) + (n3 - 1))[::-1]
-                wheel_node_idx[0, n1 + n2] = start_num * single_node_num + 2 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
-                wheel_node_idx[0, n1 + n2 + 1:n1 + n2 + n3] = np.arange(
+                # wheel_node_idx[0, 0] = start_num * single_node_num + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
+                # wheel_node_idx[0, 1:n1] = np.arange(
+                #     start_num * single_node_num + number_of_key_points + (n1 - 1) + (n2 - 1) + (n3 - 1),
+                #     start_num * single_node_num + number_of_key_points + 2 * (n1 - 1) + (n2 - 1) + (n3 - 1))[::-1]
+                # wheel_node_idx[0, n1] = start_num * single_node_num + 1 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
+                # wheel_node_idx[0, n1 + 1:n1 + n2] = np.arange(
+                #     start_num * single_node_num + number_of_key_points + 2 * (n1 - 1) + (n2 - 1) + (n3 - 1),
+                #     start_num * single_node_num + number_of_key_points + 2 * (n1 - 1) + 2 * (n2 - 1) + (n3 - 1))[::-1]
+                #
+                wheel_node_idx[0, 0] = start_num * single_node_num + 2 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
+                wheel_node_idx[0, 1:n3] = np.arange(
                     start_num * single_node_num + number_of_key_points + 2 * (n1 - 1) + 2 * (n2 - 1) + (n3 - 1),
                     start_num * single_node_num + number_of_key_points + 2 * (n1 - 1) + 2 * (n2 - 1) + 2 * (n3 - 1))[::-1]
-                wheel_node_idx[0, n1 + n2 + n3] = start_num * single_node_num + 3 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
+                wheel_node_idx[0, n3] = start_num * single_node_num + 3 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
                 if target_tooth_tag[0] == 0 and (target_tooth_tag[-1] != (z - 1)):
-                    wheel_node_idx[-1, 0] = 0
-                    wheel_node_idx[-1, 1:n1] = np.arange(number_of_key_points, number_of_key_points + n1 - 1)
-                    wheel_node_idx[-1, n1] = 1
-                    wheel_node_idx[-1, n1 + 1:n1 + n2] = np.arange(number_of_key_points + (n1 - 1),
-                                                                    number_of_key_points + (n1 - 1) + (n2 - 1))
-                    wheel_node_idx[-1, n1 + n2] = 2
-                    wheel_node_idx[-1, n1 + n2 + 1:n1 + n2 + n3] = np.arange(number_of_key_points + (n1 - 1) + (n2 - 1),
+                    # wheel_node_idx[-1, 0] = 0
+                    # wheel_node_idx[-1, 1:n1] = np.arange(number_of_key_points, number_of_key_points + n1 - 1)
+                    # wheel_node_idx[-1, n1] = 1
+                    # wheel_node_idx[-1, n1 + 1:n1 + n2] = np.arange(number_of_key_points + (n1 - 1),
+                    #                                                 number_of_key_points + (n1 - 1) + (n2 - 1))
+                    wheel_node_idx[-1, 0] = 2
+                    wheel_node_idx[-1, 1:n3] = np.arange(number_of_key_points + (n1 - 1) + (n2 - 1),
                                                                    number_of_key_points + (n1 - 1) + (n2 - 1) + (n3 - 1))
-                    wheel_node_idx[-1, n1 + n2 + n3] = 3
+                    wheel_node_idx[-1, n3] = 3
                 elif target_tooth_tag[0] != 0:
-                    wheel_node_idx[-1, 0] = 0
-                    wheel_node_idx[-1, 1:n1] = np.arange(4, 4 + n1 - 1)[::-1]
-                    wheel_node_idx[-1, n1] = 1
-                    wheel_node_idx[-1, n1 + 1:n1 + n2] = np.arange(4 + (n1 - 1), 4 + (n1 - 1) + (n2 - 1))[::-1]
-                    wheel_node_idx[-1, n1 + n2] = 2
-                    wheel_node_idx[-1, n1 + n2 + 1:n1 + n2 + n3] = np.arange(4 + (n1 - 1) + (n2 - 1), 4 + (n1 - 1) + (n2 - 1) + (n3-1))[::-1]
-                    wheel_node_idx[-1, n1 + n2 + n3] = 3
+                    # wheel_node_idx[-1, 0] = 0
+                    # wheel_node_idx[-1, 1:n1] = np.arange(4, 4 + n1 - 1)[::-1]
+                    # wheel_node_idx[-1, n1] = 1
+                    # wheel_node_idx[-1, n1 + 1:n1 + n2] = np.arange(4 + (n1 - 1), 4 + (n1 - 1) + (n2 - 1))[::-1]
+                    wheel_node_idx[-1, 0] = 2
+                    wheel_node_idx[-1, 1:n3] = np.arange(4 + (n1 - 1) + (n2 - 1), 4 + (n1 - 1) + (n2 - 1) + (n3-1))[::-1]
+                    wheel_node_idx[-1, n3] = 3
                 else:
-                    wheel_node_idx[-1, 0] = end_num * single_node_num + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
-                    wheel_node_idx[-1, 1:n1] = np.arange(
-                        end_num * single_node_num + 8 + (n1 - 1) + (n2 - 1) + (n3 - 1),
-                        end_num * single_node_num + 8 + 2*(n1 - 1) + (n2 - 1) + (n3 - 1))[::-1]
-                    wheel_node_idx[-1, n1] = end_num * single_node_num + 1 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
-                    wheel_node_idx[-1, n1 + 1:n1 + n2] = np.arange(
-                        end_num * single_node_num + 8 + 2*(n1 - 1) + (n2 - 1) + (n3 - 1),
-                        end_num * single_node_num + 8 + 2*(n1 - 1) + 2*(n2 - 1) + (n3 - 1))[::-1]
-                    wheel_node_idx[-1, n1 + n2] = end_num * single_node_num + 2 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
-                    wheel_node_idx[-1, n1 + n2 + 1:n1 + n2 + n3] = np.arange(
+                    # wheel_node_idx[-1, 0] = end_num * single_node_num + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
+                    # wheel_node_idx[-1, 1:n1] = np.arange(
+                    #     end_num * single_node_num + 8 + (n1 - 1) + (n2 - 1) + (n3 - 1),
+                    #     end_num * single_node_num + 8 + 2*(n1 - 1) + (n2 - 1) + (n3 - 1))[::-1]
+                    # wheel_node_idx[-1, n1] = end_num * single_node_num + 1 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
+                    # wheel_node_idx[-1, n1 + 1:n1 + n2] = np.arange(
+                    #     end_num * single_node_num + 8 + 2*(n1 - 1) + (n2 - 1) + (n3 - 1),
+                    #     end_num * single_node_num + 8 + 2*(n1 - 1) + 2*(n2 - 1) + (n3 - 1))[::-1]
+                    wheel_node_idx[-1, 0] = end_num * single_node_num + 2 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
+                    wheel_node_idx[-1, 1:n3] = np.arange(
                         end_num * single_node_num + 8 + 2*(n1 - 1) + 2*(n2 - 1) + (n3 - 1),
                         end_num * single_node_num + 8 + 2 * (n1 - 1) + 2 * (n2 - 1) + 2*(n3 - 1))[::-1]
-                    wheel_node_idx[-1, n1 + n2 + n3] = end_num * single_node_num + 3 + 4 + (n1 - 1) + (n2 - 1) + (
-                                n3 - 1)
+                    wheel_node_idx[-1, n3] = end_num * single_node_num + 3 + 4 + (n1 - 1) + (n2 - 1) + (n3 - 1)
 
                 wheel_cell = np.zeros((len(phi) + 1, len(r) - 1, 4), dtype=np.int64)
 

diff --git a/app/gearx/test/test_gear_system.py b/app/gearx/test/test_gear_system.py
@@ -35,7 +35,7 @@ def test_external_gear(self):
         tooth_width = data['tooth_width']
         inner_diam = data['inner_diam']  # 轮缘内径
         # m_n = 2.25
-        z = 18
+        # z = 18
         # alpha_n = 17.5
         # beta = 30
         # x_n = 0.4
@@ -424,44 +424,11 @@ def test_get_part_internal_gear(self):
 
         hex_mesh = internal_gear.generate_hexahedron_mesh()
         target_hex_mesh = internal_gear.set_target_tooth([0, 1, 2, 78, 79])
-        # part_external_gear_mesh = internal_gear.set_target_tooth([0, 1, 2, 78, 79])
-        # part_external_gear_mesh_wheel = internal_gear.set_target_tooth([0, 1, 2, 78, 79], get_wheel=True)
-        #
-        # part_external_gear_mesh.to_vtk(fname='../data/part_internal_gear_mesh.vtu')
-        # part_external_gear_mesh_wheel.to_vtk(fname='../data/part_internal_gear_hex_mesh_wheel.vtu')
-
-        n = 15
-        helix_d = np.linspace(internal_gear.d, internal_gear.d_a, n)
-        helix_width = np.linspace(0, internal_gear.tooth_width, n)
-        helix_node = internal_gear.cylindrical_to_cartesian(helix_d, helix_width)
-
-        target_cell_idx = np.zeros(n, np.int32)
-        face_normal = np.zeros((n, 3), np.float64)
-        parameters = np.zeros((n, 3), np.float64)
-        for i, t_node in enumerate(helix_node):
-            target_cell_idx[i], face_normal[i], parameters[i] = internal_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
-
-        # 寻找外圈上节点
-        node = target_hex_mesh.node
-        node_r = np.sqrt(node[:, 0] ** 2 + node[:, 1] ** 2)
-        is_outer_node = np.abs(node_r - internal_gear.outer_diam / 2) < 1e-11
-        outer_node_idx = np.where(np.abs(node_r - internal_gear.outer_diam / 2)<1e-11)[0]
+        part_internal_gear_mesh = internal_gear.set_target_tooth([0, 1, 2, 78, 79])
+        part_internal_gear_mesh_wheel = internal_gear.set_target_tooth([0, 1, 2, 78, 79], get_wheel=True)
 
-        with open('../data/part_internal_gear_data.pkl', 'wb') as f:
-            pickle.dump({'internal_gear': internal_gear, 'hex_mesh': target_hex_mesh, 'helix_node': helix_node,
-                         'target_cell_idx': target_cell_idx, 'parameters': parameters,
-                         'is_outer_node': is_outer_node, 'inner_node_idx': outer_node_idx}, f)
+        part_internal_gear_mesh.to_vtk(fname='../data/part_internal_gear_mesh.vtu')
+        part_internal_gear_mesh_wheel.to_vtk(fname='../data/part_internal_gear_hex_mesh_wheel_new.vtu')
 
     def test_get_part_external_gear(self):
         with open('../data/external_gear_data.json', 'r') as file:
@@ -501,8 +468,8 @@ def test_get_profile_face_normal_external(self):
         hex_mesh = data['hex_mesh']
         quad_mesh = data['quad_mesh']
 
-        mesh = external_gear.set_target_tooth([0, 1, 17])
-        external_gear.target_quad_mesh.to_vtk(fname='../data/face_normal_target_quad_mesh.vtu')
+        mesh = external_gear.set_target_tooth([0, 1, 18])
+        # external_gear.target_quad_mesh.to_vtk(fname='../data/face_normal_target_quad_mesh.vtu')
 
         GD = mesh.geo_dimension()
         NC = mesh.number_of_cells()
@@ -527,7 +494,7 @@ def test_get_profile_face_normal_internal(self):
         quad_mesh = data['quad_mesh']
 
         mesh = internal_gear.set_target_tooth([78, 79, 0, 1, 2])
-        internal_gear.target_quad_mesh.to_vtk(fname='../data/face_normal_target_quad_mesh_internal.vtu')
+        # internal_gear.target_quad_mesh.to_vtk(fname='../data/face_normal_target_quad_mesh_internal.vtu')
 
         # 齿面上的节点索引和坐标
         node_indices_tuple, noe_coord_tuple, profile_node_normal = internal_gear.get_profile_node_index(tooth_tag=0)

diff --git a/app/gearx/utils.py b/app/gearx/utils.py
@@ -722,7 +722,7 @@ def face_normal_bilinear(cell_nodes, u, v, w, error=1e-3, is_normalize=True):
 def get_face_normal_with_reference(cell_node, reference_normal, error=1e-3, is_normalize=True):
     '''
     给定一组六面体单元，以及一个参考法向，计算该单元的六个面的法向，并返回与参考法向最接近的法向
-    :param cell: 计算单元节点坐标
+    :param cell_node: 计算单元节点坐标
     :param reference_normal: 参考法向
     :param error: 误差限制
     :param is_normalize: 是否归一化

diff --git a/app/soptx/soptx/opt/mma.py b/app/soptx/soptx/opt/mma.py
@@ -282,8 +282,8 @@ def optimize(self, rho: TensorLike, **kwargs) -> Tuple[TensorLike, OptimizationH
 
             # MMA 方法
             volfrac = self.constraint.volume_fraction
-            # 体积约束对设计变量的标准化梯度
             fval = bm.sum(rho_phys[:]) / (volfrac * obj_grad.shape[0]) - 1
+            # 体积约束对设计变量的标准化梯度
             dfdx = con_grad[:, None].T / (volfrac * con_grad.shape[0]) # (m, n)
             rho_new, low, upp = self._solve_subproblem(
                                         xval=rho[:, None], fval=fval, 

diff --git a/app/soptx/soptx/pde/cantilever_2d.py b/app/soptx/soptx/pde/cantilever_2d.py
@@ -7,6 +7,9 @@
 from builtins import list
 
 class Cantilever2dData1:
+    '''
+    模型来源论文: Efficient topology optimization in MATLAB using 88 lines of code
+    '''
     def __init__(self, 
                 xmin: float, xmax: float, 
                 ymin: float, ymax: float):

diff --git a/app/soptx/soptx/pde/cantilever_3d.py b/app/soptx/soptx/pde/cantilever_3d.py
@@ -7,6 +7,9 @@
 
 
 class Cantilever3dData1:
+    '''
+    模型来源论文: An efficient 3D topology optimization code written in Matlab
+    '''
     def __init__(self,
                 xmin: float=0, xmax: float=60, 
                 ymin: float=0, ymax: float=20,

diff --git a/app/soptx/soptx/pde/mbb_beam_2d.py b/app/soptx/soptx/pde/mbb_beam_2d.py
@@ -7,6 +7,9 @@
 from builtins import list
 
 class MBBBeam2dData1:
+    '''
+    模型来源论文: Efficient topology optimization in MATLAB using 88 lines of code
+    '''
     def __init__(self, 
                 xmin: float=0, xmax: float=60, 
                 ymin: float=0, ymax: float=20):

diff --git a/app/soptx/soptx/tests/test_cantilever_2d.py b/app/soptx/soptx/tests/test_cantilever_2d.py
@@ -192,30 +192,41 @@ def run_optimization_test(config: TestConfig) -> Dict[str, Any]:
     }
 
 if __name__ == "__main__":
+    base_dir = '/home/heliang/FEALPy_Development/fealpy/app/soptx/soptx/vtu'
+
     # 使用 OC 优化器的配置
+    '''
+    参数来源论文: Efficient topology optimization in MATLAB using 88 lines of code
+    '''
+    filter_type = 'sensitivity'
+    optimizer_type = 'oc'
     config1 = TestConfig(
         nx=160, ny=100,
         volume_fraction=0.4,
         filter_radius=6.0,
-        filter_type='sensitivity',
-        max_iterations=100,
-        save_dir='/home/heliang/FEALPy_Development/fealpy/app/soptx/soptx/tests/cantilever_2d_oc',
+        filter_type=filter_type,       # 指定使用灵敏度滤波器
+        save_dir=f'{base_dir}/cantilever_2d_{filter_type}_{optimizer_type}',
         mesh_type='uniform_mesh_2d',
         assembly_method=AssemblyMethod.FAST_STRESS_UNIFORM,
-        optimizer_type='oc'  # 指定使用 OC 优化器
+        optimizer_type=optimizer_type,  # 指定使用 OC 优化器
+        max_iterations=200,
+        tolerance=0.01
     )
 
     # 使用 MMA 优化器的配置
+    filter_type = 'sensitivity'
+    optimizer_type = 'mma'
     config2 = TestConfig(
         nx=160, ny=100,
         volume_fraction=0.4,
         filter_radius=6.0,
-        filter_type='density',
-        max_iterations=100,
-        save_dir='/home/heliang/FEALPy_Development/fealpy/app/soptx/soptx/tests/cantilever_2d_mma',
+        filter_type=filter_type,       # 指定使用灵敏度滤波器
+        save_dir=f'{base_dir}/cantilever_2d_{filter_type}_{optimizer_type}',
         mesh_type='uniform_mesh_2d',
         assembly_method=AssemblyMethod.FAST_STRESS_UNIFORM,
-        optimizer_type='mma'  # 指定使用 MMA 优化器
+        optimizer_type=optimizer_type,  # 指定使用 OC 优化器
+        max_iterations=200,
+        tolerance=0.01
     )
 
     result = run_optimization_test(config2)

diff --git a/app/soptx/soptx/tests/test_cantilever_3d.py b/app/soptx/soptx/tests/test_cantilever_3d.py
@@ -192,6 +192,9 @@ def run_optimization_test(config: TestConfig) -> Dict[str, Any]:
     base_dir = '/home/heliang/FEALPy_Development/fealpy/app/soptx/soptx/vtu'
 
     # 使用 OC 优化器的配置
+    '''
+    参数来源论文: An efficient 3D topology optimization code written in Matlab
+    '''
     filter_type = 'density'
     optimizer_type = 'oc'
     config1 = TestConfig(
@@ -209,6 +212,9 @@ def run_optimization_test(config: TestConfig) -> Dict[str, Any]:
     )
 
     # 使用 MMA 优化器的配置
+    '''
+    参数来源论文: An efficient 3D topology optimization code written in Matlab
+    '''
     filter_type = 'density'
     optimizer_type = 'mma'
     config2 = TestConfig(