diff --git a/app/soptx/linear_elasticity/jy_box_api.py b/app/soptx/linear_elasticity/jy_box_api.py index 99f4dc815..a46b82dc5 100644 --- a/app/soptx/linear_elasticity/jy_box_api.py +++ b/app/soptx/linear_elasticity/jy_box_api.py @@ -214,6 +214,106 @@ def compute_element_equivalent_strain(strain, nu): space = LagrangeFESpace(mesh, p=1, ctype='C') cell2dof = space.cell_to_dof() +###########################333 +# bcs1 = (bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64),) +# p1 = mesh.bc_to_point(bc=bcs1) +# bcs2 = (bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64),) +# p2 = mesh.bc_to_point(bc=bcs2) +# bcs3 = (bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64),) +# p3 = mesh.bc_to_point(bc=bcs3) +# bcs4 = (bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64),) +# p4 = mesh.bc_to_point(bc=bcs4) +# bcs5 = (bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64),) +# p5 = mesh.bc_to_point(bc=bcs5) +# bcs6 = (bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64),) +# p6 = mesh.bc_to_point(bc=bcs6) +# bcs7 = (bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64),) +# p7 = mesh.bc_to_point(bc=bcs7) +# bcs8 = (bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64),) +# p8 = mesh.bc_to_point(bc=bcs8) + +bcs1 = (bm.tensor([[1, 0]], dtype=bm.float64), + bm.tensor([[1, 0]], dtype=bm.float64), + bm.tensor([[1, 0]], dtype=bm.float64),) +p1 = mesh.bc_to_point(bc=bcs1) +bcs5 = (bm.tensor([[0, 1]], dtype=bm.float64), + bm.tensor([[1, 0]], dtype=bm.float64), + bm.tensor([[1, 0]], dtype=bm.float64),) +p5 = mesh.bc_to_point(bc=bcs5) +bcs7 = (bm.tensor([[0, 1]], dtype=bm.float64), + bm.tensor([[0, 1]], dtype=bm.float64), + bm.tensor([[1, 0]], dtype=bm.float64),) +p7 = mesh.bc_to_point(bc=bcs7) +bcs3 = (bm.tensor([[1, 0]], dtype=bm.float64), + bm.tensor([[0, 1]], dtype=bm.float64), + bm.tensor([[1, 0]], dtype=bm.float64),) +p3 = mesh.bc_to_point(bc=bcs3) +bcs2 = (bm.tensor([[1, 0]], dtype=bm.float64), + bm.tensor([[1, 0]], dtype=bm.float64), + bm.tensor([[0, 1]], dtype=bm.float64),) +p2 = mesh.bc_to_point(bc=bcs2) +bcs6 = (bm.tensor([[0, 1]], dtype=bm.float64), + bm.tensor([[1, 0]], dtype=bm.float64), + bm.tensor([[0, 1]], dtype=bm.float64),) +p6 = mesh.bc_to_point(bc=bcs6) +bcs8 = (bm.tensor([[0, 1]], dtype=bm.float64), + bm.tensor([[0, 1]], dtype=bm.float64), + bm.tensor([[0, 1]], dtype=bm.float64),) +p8 = mesh.bc_to_point(bc=bcs8) +bcs4 = (bm.tensor([[1, 0]], dtype=bm.float64), + bm.tensor([[0, 1]], dtype=bm.float64), + bm.tensor([[0, 1]], dtype=bm.float64),) +p4 = mesh.bc_to_point(bc=bcs4) + +test0 = space.function() +test0[0] = 1 +print(f"1:{p1}, {test0(bcs1)}, {test0(bcs2)}, {test0(bcs3)}, {test0(bcs4)}, {test0(bcs5)}, {test0(bcs6)}, {test0(bcs7)}, {test0(bcs8)}") + +test1 = space.function() +test1[1] = 1 +print(f"2:{p2}, {test1(bcs2)}, {test1(bcs1)}, {test1(bcs3)}, {test1(bcs4)}, {test1(bcs5)}, {test1(bcs6)}, {test1(bcs7)}, {test1(bcs8)}") + +test2 = space.function() +test2[2] = 1 +print(f"3:{p3}, {test2(bcs3)}, {test2(bcs1)}, {test2(bcs2)}, {test2(bcs4)}, {test2(bcs5)}, {test2(bcs6)}, {test2(bcs7)}, {test2(bcs8)}") + +test3 = space.function() +test3[3] = 1 +print(f"4:{p4}, {test3(bcs4)}, {test3(bcs1)}, {test3(bcs2)}, {test3(bcs3)}, {test3(bcs5)}, {test3(bcs6)}, {test3(bcs7)}, {test3(bcs8)}") + +test4 = space.function() +test4[4] = 1 +print(f"5:{p5}, {test4(bcs5)}, {test4(bcs1)}, {test4(bcs2)}, {test4(bcs3)}, {test4(bcs4)}, {test4(bcs6)}, {test4(bcs7)}, {test4(bcs8)}") + +test5 = space.function() +test5[5] = 1 +print(f"6:{p6}, {test5(bcs6)}, {test5(bcs1)}, {test5(bcs2)}, {test5(bcs3)}, {test5(bcs4)}, {test5(bcs5)}, {test5(bcs7)}, {test5(bcs8)}") + +test6 = space.function() +test6[6] = 1 +print(f"7:{p7}, {test6(bcs7)}, {test6(bcs1)}, {test6(bcs2)}, {test6(bcs3)}, {test6(bcs4)}, {test6(bcs5)}, {test6(bcs6)}, {test6(bcs8)}") + +test7 = space.function() +test7[7] = 1 +print(f"8:{p8}, {test7(bcs8)}, {test7(bcs1)}, {test7(bcs2)}, {test7(bcs3)}, {test7(bcs4)}, {test7(bcs5)}, {test7(bcs6)}, {test7(bcs7)}") +##############3 + q = 2 # tensor_space = TensorFunctionSpace(space, shape=(-1, 3)) # gd_priority diff --git a/app/soptx/linear_elasticity/jy_gear_api.py b/app/soptx/linear_elasticity/jy_gear_api.py index eea6188fa..34bb758b4 100644 --- a/app/soptx/linear_elasticity/jy_gear_api.py +++ b/app/soptx/linear_elasticity/jy_gear_api.py @@ -38,6 +38,25 @@ def compute_equivalent_strain(strain, nu): return equiv_strain +def compute_equivalent_stress(stress, nu): + sxx = stress[..., 0, 0] + syy = stress[..., 1, 1] + szz = stress[..., 2, 2] + sxy = stress[..., 0, 1] + syz = stress[..., 1, 2] + sxz = stress[..., 0, 2] + + d1 = sxx - syy + d2 = syy - szz + d3 = szz - sxx + + # Combine all terms + equiv_stress = (d1**2 + d2**2 + d3**2 + 6.0 * (sxy**2 + syz**2 + sxz**2)) + + equiv_stress = bm.sqrt(equiv_stress / 2.0) + + return equiv_stress + with open('/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/external_gear_data_part.pkl', 'rb') as f: data = pickle.load(f) @@ -214,6 +233,31 @@ def compute_equivalent_strain(strain, nu): for c in range(NC): uh_cell[c] = uh[cell2tdof[c]] +# bcs1 = (bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64),) +# bcs5 = (bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64),) +# bcs7 = (bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64),) +# bcs3 = (bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64),) +# bcs2 = (bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64),) +# bcs6 = (bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64),) +# bcs8 = (bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64),) +# bcs4 = (bm.tensor([[1, 0]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64), +# bm.tensor([[0, 1]], dtype=bm.float64),) + bcs1 = (bm.tensor([[1, 0]], dtype=bm.float64), bm.tensor([[1, 0]], dtype=bm.float64), bm.tensor([[1, 0]], dtype=bm.float64),) @@ -256,18 +300,36 @@ def compute_equivalent_strain(strain, nu): strain_xx = strain[..., 0, 0] # (NC, 8) strain_yy = strain[..., 1, 1] # (NC, 8) strain_zz = strain[..., 2, 2] # (NC, 8) -strain_xy = strain[..., 0, 1] # (NC, 8) +strain_xy = strain[..., 0, 1] # (NC, 8) strain_yz = strain[..., 1, 2] # (NC, 8) strain_xz = strain[..., 0, 2] # (NC, 8) +trace_e = bm.einsum("...ii", strain) # (NC, 8) +I = bm.eye(GD, dtype=bm.float64) +stress = 2 * mu * strain + lam * trace_e[..., None, None] * I # (NC, 8, GD, GD) + +stress_xx = stress[..., 0, 0] # (NC, 8) +stress_yy = stress[..., 1, 1] # (NC, 8) +stress_zz = stress[..., 2, 2] # (NC, 8) +stress_xy = stress[..., 0, 1] # (NC, 8) +stress_yz = stress[..., 1, 2] # (NC, 8) +stress_xz = stress[..., 0, 2] # (NC, 8) + equiv_strain = compute_equivalent_strain(strain, nu) # (NC, 8) -nodal_equiv_strain = bm.zeros(NN, dtype=bm.float64) +equiv_stress = compute_equivalent_stress(stress, nu) # (NC, 8) + num_count = bm.zeros(NN, dtype=bm.int32) -bm.add_at(nodal_equiv_strain, cell2dof.flatten(), equiv_strain.flatten()) bm.add_at(num_count, cell2dof.flatten(), 1) + +nodal_equiv_strain = bm.zeros(NN, dtype=bm.float64) +bm.add_at(nodal_equiv_strain, cell2dof.flatten(), equiv_strain.flatten()) nodal_equiv_strain = nodal_equiv_strain / num_count +nodal_equiv_stress = bm.zeros(NN, dtype=bm.float64) +bm.add_at(nodal_equiv_stress, cell2dof.flatten(), equiv_stress.flatten()) +nodal_equiv_stress = nodal_equiv_stress / num_count + nodal_strain_xx = bm.zeros(NN, dtype=bm.float64) bm.set_at(nodal_strain_xx, cell2dof.flatten(), strain_xx.flatten()) nodal_strain_yy = bm.zeros(NN, dtype=bm.float64) @@ -281,14 +343,79 @@ def compute_equivalent_strain(strain, nu): nodal_strain_xz = bm.zeros(NN, dtype=bm.float64) bm.set_at(nodal_strain_xz, cell2dof.flatten(), strain_xz.flatten()) +nodal_stress_xx = bm.zeros(NN, dtype=bm.float64) +bm.set_at(nodal_stress_xx, cell2dof.flatten(), stress_xx.flatten()) +nodal_stress_yy = bm.zeros(NN, dtype=bm.float64) +bm.set_at(nodal_stress_yy, cell2dof.flatten(), stress_yy.flatten()) +nodal_stress_zz = bm.zeros(NN, dtype=bm.float64) +bm.set_at(nodal_stress_zz, cell2dof.flatten(), stress_zz.flatten()) +nodal_stress_xy = bm.zeros(NN, dtype=bm.float64) +bm.set_at(nodal_stress_xy, cell2dof.flatten(), stress_xy.flatten()) +nodal_stress_yz = bm.zeros(NN, dtype=bm.float64) +bm.set_at(nodal_stress_yz, cell2dof.flatten(), stress_yz.flatten()) +nodal_stress_xz = bm.zeros(NN, dtype=bm.float64) +bm.set_at(nodal_stress_xz, cell2dof.flatten(), stress_xz.flatten()) + +mesh.celldata['equiv_strain'] = equiv_strain[:] +mesh.celldata['equiv_stress'] = equiv_stress[:] +mesh.celldata['strain_xx'] = strain_xx[:] +mesh.celldata['strain_yy'] = strain_yy[:] +mesh.celldata['strain_zz'] = strain_zz[:] +mesh.celldata['strain_xy'] = strain_xy[:] +mesh.celldata['strain_yz'] = strain_yz[:] +mesh.celldata['strain_xz'] = strain_xz[:] +mesh.celldata['stress_xx'] = stress_xx[:] +mesh.celldata['stress_yy'] = stress_yy[:] +mesh.celldata['stress_zz'] = stress_zz[:] +mesh.celldata['stress_xy'] = stress_xy[:] +mesh.celldata['stress_yz'] = stress_yz[:] +mesh.celldata['stress_xz'] = stress_xz[:] mesh.nodedata['nodal_equiv_strain'] = nodal_equiv_strain[:] +mesh.nodedata['nodal_equiv_stress'] = nodal_equiv_stress[:] mesh.nodedata['nodal_strain_xx'] = nodal_strain_xx[:] mesh.nodedata['nodal_strain_yy'] = nodal_strain_yy[:] mesh.nodedata['nodal_strain_zz'] = nodal_strain_zz[:] mesh.nodedata['nodal_strain_xy'] = nodal_strain_xy[:] mesh.nodedata['nodal_strain_yz'] = nodal_strain_yz[:] mesh.nodedata['nodal_strain_xz'] = nodal_strain_xz[:] +mesh.nodedata['nodal_stress_xx'] = nodal_stress_xx[:] +mesh.nodedata['nodal_stress_yy'] = nodal_stress_yy[:] +mesh.nodedata['nodal_stress_zz'] = nodal_stress_zz[:] +mesh.nodedata['nodal_stress_xy'] = nodal_stress_xy[:] +mesh.nodedata['nodal_stress_yz'] = nodal_stress_yz[:] +mesh.nodedata['nodal_stress_xz'] = nodal_stress_xz[:] mesh.to_vtk('/home/heliang/FEALPy_Development/fealpy/app/soptx/linear_elasticity/gear_fealpy.vtu') - +# # 首先将应变张量转换为Voigt记号形式 +# strain_voigt = bm.zeros((NC, 8, 6), dtype=bm.float64) +# strain_voigt[..., 0] = strain[..., 0, 0] # εxx +# strain_voigt[..., 1] = strain[..., 1, 1] # εyy +# strain_voigt[..., 2] = strain[..., 2, 2] # εzz +# strain_voigt[..., 3] = 2 * strain[..., 0, 1] # 2εxy +# strain_voigt[..., 4] = 2 * strain[..., 1, 2] # 2εyz +# strain_voigt[..., 5] = 2 * strain[..., 0, 2] # 2εxz + +# error_xx = strain_xx - strain_voigt[..., 0] +# error_yy = strain_yy - strain_voigt[..., 1] +# error_zz = strain_zz - strain_voigt[..., 2] +# error_xy = strain_xy - strain_voigt[..., 3] / 2 +# error_yz = strain_yz - strain_voigt[..., 4] / 2 +# error_xz = strain_xz - strain_voigt[..., 5] / 2 + +# # 使用弹性矩阵计算应力（Voigt记号形式） +# D = linear_elastic_material.elastic_matrix() # (1, 1, 6, 6) +# D_expanded = bm.broadcast_to(D, (NC, 8, 6, 6)) +# stress_voigt = bm.einsum('cqij, cqj -> cqi', D_expanded, strain_voigt) # (NC, 8, 6) + +# # 将 Voigt 记号形式的应力转换回张量形式 +# stress = bm.zeros((NC, 8, GD, GD), dtype=bm.float64) +# stress[..., 0, 0] = stress_voigt[..., 0] # σxx +# stress[..., 1, 1] = stress_voigt[..., 1] # σyy +# stress[..., 2, 2] = stress_voigt[..., 2] # σzz +# stress[..., 0, 1] = stress_voigt[..., 3] / 2 # σxy +# stress[..., 1, 0] = stress_voigt[..., 3] / 2 # σyx +# stress[..., 1, 2] = stress_voigt[..., 4] / 2 # σyz +# stress[..., 2, 1] = stress_voigt[..., 4] / 2 # σzy +# stress[..., 0, 2] = stress_voigt[..., 5] / 2 # σxz +# stress[..., 2, 0] = stress_voigt[..., 5] / 2 # σzx print("-----------") \ No newline at end of file