diff --git a/app/soptx/linear_elasticity/JingYiGearProject/external_gear_profile.py b/app/soptx/linear_elasticity/JingYiGearProject/external_gear_profile.py index c48c1b967..7d0faf6ec 100644 --- a/app/soptx/linear_elasticity/JingYiGearProject/external_gear_profile.py +++ b/app/soptx/linear_elasticity/JingYiGearProject/external_gear_profile.py @@ -1,5 +1,5 @@ """ -外齿轮 15 个载荷点的算例 +外齿轮 1 个载荷点的算例 (左右齿面) """ from fealpy.backend import backend_manager as bm from fealpy.mesh import HexahedronMesh @@ -14,21 +14,14 @@ from soptx.utils import timer from app.soptx.linear_elasticity.JingYiGearProject.utils import export_to_inp -from app.gearx.gear import ExternalGear, InternalGear +from app.gearx.gear import ExternalGear import json -def compute_strain_stress(tensor_space, uh, B_BBar, D): - cell2tdof = tensor_space.cell_to_dof() - cuh = uh[cell2tdof] # (NC, TLDOF) - strain = bm.einsum('cqil, cl -> cqi', B_BBar, cuh) # (NC, NQ, 6) - stress = bm.einsum('cqij, cqi -> cqj', D, strain) # (NC, NQ, 6) - - return strain, stress - bm.set_backend('numpy') -with open('/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/JingYiGearProject/json/external_gear_data.json', 'r') \ - as file:data = json.load(file) +with open( +'/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/JingYiGearProject/json/external_gear_data.json', 'r') \ + as file:data = json.load(file) m_n = data['mn'] # 法向模数 z = data['z'] # 齿数 alpha_n = data['alpha_n'] # 法向压力角 @@ -99,22 +92,23 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D): bform.add_integrator(integrator_K) K = bform.assembly(format='csr') -# 齿面上的节点索引和坐标 -node_indices_tuple, noe_coord_tuple = external_gear.get_profile_node_index(tooth_tag=0) +# 齿面上的节点索引, 坐标和内法线方向 +node_indices_tuple, node_coord_tuple, profile_node_normal_tuple = external_gear.get_profile_node_index(tooth_tag=0) node_indices_left = node_indices_tuple[0].reshape(-1, 1) node_indices_right = node_indices_tuple[1].reshape(-1, 1) node_indices = bm.concatenate([node_indices_left, node_indices_right], axis=0) # (NPN, 1) -node_coord_left = noe_coord_tuple[0].reshape(-1, 3) -node_coord_right = noe_coord_tuple[1].reshape(-1, 3) +node_coord_left = node_coord_tuple[0].reshape(-1, 3) +node_coord_right = node_coord_tuple[1].reshape(-1, 3) node_coord = bm.concatenate([node_coord_left, node_coord_right], axis=0) # (NPN, GD) +profile_node_normal_left = profile_node_normal_tuple[0].reshape(-1, 3) +profile_node_normal_right = profile_node_normal_tuple[1].reshape(-1, 3) +profile_node_normal = bm.concatenate([profile_node_normal_left, profile_node_normal_right], axis=0) # (NPN, GD) # 齿面上的节点数 NPN = node_indices.shape[0] -# TODO 齿面上节点的内法线方向 -face_normal = bm.ones((NPN, 3), bm.float64) # 节点载荷值 load_values = 1000 # 所有节点的内法线载荷向量 -P = load_values * face_normal # (NPN, 3) +P = load_values * profile_node_normal # (NPN, GD) # 齿面上节点对应的全局自由度编号（跟顺序有关） if tensor_space.dof_priority: dof_indices = bm.stack([sgdof * d + node_indices.reshape(-1) for d in range(GD)], axis=1) # (NPN, GD) @@ -134,7 +128,7 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D): method='interp') # 齿面上所有节点的位移结果 uh_profiles = bm.zeros((NPN, GD), dtype=bm.float64) # (NPN, GD) -for i in range(3): +for i in range(NPN): # 创建计时器 t = timer(f"Timing_{i}") next(t) # 启动计时器 @@ -153,10 +147,10 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D): logger.setLevel('INFO') uh = tensor_space.function() # uh[:] = cg(K, F, maxiter=10000, atol=1e-8, rtol=1e-8) - uh[:] = spsolve(K, F, solver="mumps") + uh[:] = spsolve(K, F, solver="mumps") # (tgdof, ) t.send('求解时间') # 获取齿面上节点的位移 - uh_profile = uh[dof_indices[i, :]] + uh_profile = uh[dof_indices[i, :]] # (GD, ) uh_profiles[i, :] = uh_profile # 计算残差向量和范数 residual = K.matmul(uh[:]) - F @@ -185,5 +179,5 @@ def compute_strain_stress(tensor_space, uh, B_BBar, D): young_modulus=206e3, poisson_ratio=0.3, density=7.85e-9, used_app='abaqus', mesh_type='hex') # 计算齿面上节点的内法线方向位移 -uh_normal = bm.sum(uh_profiles * face_normal, axis=1) # (NPN, ) +uh_normal = bm.sum(uh_profiles * profile_node_normal, axis=1) # (NPN, ) print("-----------") \ No newline at end of file diff --git a/app/soptx/linear_elasticity/JingYiGearProject/internal_gear_helix.py b/app/soptx/linear_elasticity/JingYiGearProject/internal_gear_helix.py index 4f3b5765d..c314b7ffb 100644 --- a/app/soptx/linear_elasticity/JingYiGearProject/internal_gear_helix.py +++ b/app/soptx/linear_elasticity/JingYiGearProject/internal_gear_helix.py @@ -1,5 +1,5 @@ """ -外齿轮 15 个载荷点的算例 +内齿轮 15 个载荷点的算例 """ from fealpy.backend import backend_manager as bm from fealpy.functionspace import LagrangeFESpace, TensorFunctionSpace diff --git a/app/soptx/linear_elasticity/JingYiGearProject/internal_gear_profile.py b/app/soptx/linear_elasticity/JingYiGearProject/internal_gear_profile.py new file mode 100644 index 000000000..551c79f2c --- /dev/null +++ b/app/soptx/linear_elasticity/JingYiGearProject/internal_gear_profile.py @@ -0,0 +1,183 @@ +""" +内齿轮 1 个载荷点的算例 (左右齿面) +""" +from fealpy.backend import backend_manager as bm +from fealpy.functionspace import LagrangeFESpace, TensorFunctionSpace +from fealpy.sparse import COOTensor +from fealpy.fem.linear_elastic_integrator import LinearElasticIntegrator +from fealpy.material.elastic_material import LinearElasticMaterial +from fealpy.fem.bilinear_form import BilinearForm +from fealpy.fem.dirichlet_bc import DirichletBC +from fealpy.solver import cg, spsolve +from fealpy.mesh import HexahedronMesh + +from soptx.utils import timer +from app.soptx.linear_elasticity.JingYiGearProject.utils import export_to_inp +from app.gearx.gear import InternalGear +import json + +bm.set_backend('numpy') + +with open( +'/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/JingYiGearProject/json/internal_gear_data.json', 'r') \ + as file:data = json.load(file) +m_n = data['mn'] # 法向模数 +z = data['z'] # 齿数 +alpha_n = data['alpha_n'] # 法向压力角 +beta = data['beta'] # 螺旋角 +x_n = data['xn'] # 法向变位系数 +hac = data['hac'] # 齿顶高系数 +cc = data['cc'] # 顶隙系数 +rcc = data['rcc'] # 刀尖圆弧半径 +jn = data['jn'] # 法向侧隙 +n1 = data['n1'] # 渐开线分段数 +n2 = data['n2'] # 过渡曲线分段数 +n3 = data['n3'] +na = data['na'] +nf = data['nf'] +nw = data['nw'] +tooth_width = data['tooth_width'] +outer_diam = data['outer_diam'] # 轮缘内径 +z_cutter = data['z_cutter'] +xn_cutter = data['xn_cutter'] + +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() +target_hex_mesh = internal_gear.set_target_tooth([0, 1, 2, 78, 79]) + +hex_cell = target_hex_mesh.cell +hex_node = target_hex_mesh.node +# 寻找外圈上节点 +node_r = bm.sqrt(hex_node[:, 0] ** 2 + hex_node[:, 1] ** 2) +is_outer_node = bm.abs(node_r - internal_gear.outer_diam / 2) < 1e-11 +outer_node_idx = bm.where(bm.abs(node_r - internal_gear.outer_diam / 2)<1e-11)[0] +mesh = HexahedronMesh(hex_node, hex_cell) + +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') + +# 创建有限元空间 +p = 1 +q = 2 +space = LagrangeFESpace(mesh, p=p, ctype='C') +sgdof = space.number_of_global_dofs() +print(f"sgdof: {sgdof}") +cell2dof = space.cell_to_dof() +tensor_space = TensorFunctionSpace(space, shape=(-1, 3)) # gd_priority +cell2tdof = tensor_space.cell_to_dof() +tgdof = tensor_space.number_of_global_dofs() +print(f"tgdof: {tgdof}") +tldof = tensor_space.number_of_local_dofs() + +# 组装刚度矩阵 +E = 206e3 +nu = 0.3 +lam = (E * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu)) +mu = E / (2.0 * (1.0 + nu)) +linear_elastic_material = LinearElasticMaterial(name='E_nu', + elastic_modulus=E, poisson_ratio=nu, + hypo='3D', device=bm.get_device(mesh)) +# B-Bar 修正的刚度矩阵 +integrator_K = LinearElasticIntegrator(material=linear_elastic_material, + q=q, method='C3D8_BBar') +integrator_K.keep_data(True) +_, _, D, B_BBar = integrator_K.fetch_c3d8_bbar_assembly(tensor_space) +bform = BilinearForm(tensor_space) +bform.add_integrator(integrator_K) +K = bform.assembly(format='csr') + +# 齿面上的节点索引, 坐标和内法线方向 +node_indices_tuple, node_coord_tuple, profile_node_normal_tuple = internal_gear.get_profile_node_index(tooth_tag=0) +node_indices_left = node_indices_tuple[0].reshape(-1, 1) +node_indices_right = node_indices_tuple[1].reshape(-1, 1) +node_indices = bm.concatenate([node_indices_left, node_indices_right], axis=0) # (NPN, 1) +node_coord_left = node_coord_tuple[0].reshape(-1, 3) +node_coord_right = node_coord_tuple[1].reshape(-1, 3) +node_coord = bm.concatenate([node_coord_left, node_coord_right], axis=0) # (NPN, GD) +profile_node_normal_left = profile_node_normal_tuple[0].reshape(-1, 3) +profile_node_normal_right = profile_node_normal_tuple[1].reshape(-1, 3) +profile_node_normal = bm.concatenate([profile_node_normal_left, profile_node_normal_right], axis=0) # (NPN, GD) +# 齿面上的节点数 +NPN = node_indices.shape[0] +# 节点载荷值 +load_values = 1000 +# 所有节点的内法线载荷向量 +P = load_values * profile_node_normal # (NPN, GD) +# 齿面上节点对应的全局自由度编号（跟顺序有关） +if tensor_space.dof_priority: + dof_indices = bm.stack([sgdof * d + node_indices.reshape(-1) for d in range(GD)], axis=1) # (NPN, GD) +else: + dof_indices = bm.stack([node_indices.reshape(-1) * GD + d for d in range(GD)], axis=1) # (NPN, GD) +# inp 文件中需要的固定节点索引 +fixed_node_index = bm.where(is_outer_node)[0] +# Dirichlet 边界条件 +scalar_is_bd_dof = bm.zeros(sgdof, dtype=bm.bool) +scalar_is_bd_dof[:NN] = is_outer_node +tensor_is_bd_dof = tensor_space.is_boundary_dof( + threshold=(scalar_is_bd_dof, scalar_is_bd_dof, scalar_is_bd_dof), + method='interp') +dbc = DirichletBC(space=tensor_space, + gd=bm.zeros(tgdof), + threshold=tensor_is_bd_dof, + method='interp') +# 齿面上所有节点的位移结果 +uh_profiles = bm.zeros((NPN, GD), dtype=bm.float64) # (NPN, GD) +for i in range(NPN): + # 创建计时器 + t = timer(f"Timing_{i}") + next(t) # 启动计时器 + PP = P[i, :] + # 全局载荷向量 + F = COOTensor(indices = bm.empty((1, 0), dtype=bm.int32, device=bm.get_device(space)), + values = bm.empty((0, ), dtype=bm.float64, device=bm.get_device(space)), + spshape = (tgdof, )) + indices = dof_indices[i].reshape(1, -1) + F = F.add(COOTensor(indices, PP.reshape(-1), (tgdof, ))).to_dense() # (tgdof, ) + # 处理 Dirichlet 边界条件 + K, F = dbc.apply(K, F) + t.send('准备时间') + # 计算位移 + from fealpy import logger + logger.setLevel('INFO') + uh = tensor_space.function() + # uh[:] = cg(K, F, maxiter=10000, atol=1e-8, rtol=1e-8) + uh[:] = spsolve(K, F, solver="mumps") # (tgdof, ) + t.send('求解时间') + # 获取齿面上节点的位移 + uh_profile = uh[dof_indices[i, :]] # (GD, ) + uh_profiles[i, :] = uh_profile + # 计算残差向量和范数 + residual = K.matmul(uh[:]) - F + residual_norm = bm.sqrt(bm.sum(residual * residual)) + print(f"Final residual norm: {residual_norm:.6e}") + t.send('后处理时间') + t.send(None) + if i == 0: + # 保存单个节点的位移结果 + if tensor_space.dof_priority: + uh_show = uh.reshape(GD, NN).T + else: + uh_show = uh.reshape(NN, GD) + uh_magnitude = bm.linalg.norm(uh_show, axis=1) + mesh.nodedata['uh'] = uh_show[:] + mesh.nodedata['uh_magnitude'] = uh_magnitude[:] + mesh.to_vtk(f'/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/JingYiGearProject/vtu/internal_gear_profile_fealpy_{i}.vtu') + + # 单个节点载荷的索引 + load_node_indices = node_indices[i].reshape(-1) # (1, ) + # 从全局载荷向量中提取单个载荷节点处的值 + F_load_nodes = F[dof_indices[i, :].reshape(1, -1)] # (1, GD) + export_to_inp(filename=f'/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/JingYiGearProject/inp/internal_gear_profile_abaqus_{i}.inp', + nodes=node, elements=cell, + fixed_nodes=fixed_node_index, load_nodes=load_node_indices, loads=F_load_nodes, + young_modulus=206e3, poisson_ratio=0.3, density=7.85e-9, + used_app='abaqus', mesh_type='hex') +# 计算齿面上节点的内法线方向位移 +uh_normal = bm.sum(uh_profiles * profile_node_normal, axis=1) # (NPN, ) +print("-----------") \ No newline at end of file diff --git a/app/soptx/soptx/opt/mma.py b/app/soptx/soptx/opt/mma.py index eae134c03..51d3971b9 100644 --- a/app/soptx/soptx/opt/mma.py +++ b/app/soptx/soptx/opt/mma.py @@ -120,9 +120,9 @@ def _update_asymptotes(self, else: factor = bm.ones((xval.shape[0], 1)) xxx = (xval - xold1) * (xold1 - xold2) - factor[xxx > 0] = asyincr - factor[xxx < 0] = asydecr - factor[xxx == 0] = 1.0 + epsilon = 1e-12 + factor[xxx > epsilon] = asyincr + factor[xxx < -epsilon] = asydecr self._low = xval - factor * (xold1 - self._low) self._upp = xval + factor * (self._upp - xold1) @@ -136,7 +136,7 @@ def _update_asymptotes(self, self._low = bm.minimum(self._low, lowmax) self._upp = bm.minimum(self._upp, uppmax) self._upp = bm.maximum(self._upp, uppmin) - + return self._low, self._upp def _solve_subproblem(self, @@ -151,6 +151,7 @@ def _solve_subproblem(self, """求解 MMA 子问题 Paramters + - xval (n, 1): 当前设计变量 - df0dx (n, 1): 目标函数的梯度 - dfdx (m, n): 约束函数的梯度 @@ -282,14 +283,15 @@ 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( - rho[:, None], fval=con_val, + xval=rho[:, None], fval=fval, df0dx=obj_grad[:, None], dfdx=dfdx, low=low, upp=upp, xold1=xold1[:, None], xold2=xold2[..., None] ) - + # 更新物理密度 if self.filter is not None: rho_phys = self.filter.filter_density(rho_new, filter_params) diff --git a/app/soptx/soptx/tests/test_cantilever_3d.py b/app/soptx/soptx/tests/test_cantilever_3d.py index f3e9abf01..21f939512 100644 --- a/app/soptx/soptx/tests/test_cantilever_3d.py +++ b/app/soptx/soptx/tests/test_cantilever_3d.py @@ -46,7 +46,7 @@ class TestConfig: # 新增优化器类型参数 optimizer_type: Literal['oc', 'mma'] = 'oc' # 默认使用 OC 方法 max_iterations: int = 200 - tolerance: float = 0.001 + tolerance: float = 0.01 # OC 优化器的参数 move_limit: float = 0.2 @@ -189,30 +189,39 @@ def run_optimization_test(config: TestConfig) -> Dict[str, Any]: } if __name__ == "__main__": + base_dir = '/home/heliang/FEALPy_Development/fealpy/app/soptx/soptx/vtu' + # 使用 OC 优化器的配置 + filter_type = 'density' + optimizer_type = 'oc' config1 = TestConfig( - nx=160, ny=100, nz=4, - 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_3d_oc', - mesh_type='uniform_mesh_2d', - assembly_method=AssemblyMethod.FAST_STRESS_UNIFORM, - optimizer_type='oc' # 指定使用 OC 优化器 + nx=60, ny=20, nz=4, + volume_fraction=0.3, + filter_radius=1.5, + filter_type=filter_type, # 指定使用密度滤波器 + save_dir=f'{base_dir}/cantilever_3d_{filter_type}_{optimizer_type}', + mesh_type='uniform_mesh_3d', + assembly_method=AssemblyMethod.FAST_3D_UNIFORM, + optimizer_type=optimizer_type, # 指定使用 OC 优化器 + max_iterations=200, + tolerance=0.01 + ) # 使用 MMA 优化器的配置 + filter_type = 'density' + optimizer_type = 'mma' config2 = TestConfig( nx=60, ny=20, nz=4, volume_fraction=0.3, filter_radius=1.5, - filter_type='density', - max_iterations=15, - save_dir='/home/heliang/FEALPy_Development/fealpy/app/soptx/soptx/tests/cantilever_3d_mma', + filter_type=filter_type, # 指定使用密度滤波器 + save_dir=f'{base_dir}/cantilever_3d_{filter_type}_{optimizer_type}', mesh_type='uniform_mesh_3d', assembly_method=AssemblyMethod.FAST_3D_UNIFORM, - optimizer_type='mma' # 指定使用 MMA 优化器 + optimizer_type=optimizer_type, # 指定使用 OC 优化器 + max_iterations=200, + tolerance=0.01 ) result = run_optimization_test(config2)