weihuayi · AlbertZyy · Dec 6, 2024 · Dec 5, 2024 · Dec 5, 2024 · Dec 5, 2024
diff --git a/.DS_Store b/.DS_Store
diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml
@@ -1,30 +1,36 @@
-# This workflows will upload a Python Package using Twine when a release is created
-# For more information see: https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#publishing-to-package-registries
-name: Upload Python Package
+name: Publish to PyPI
 
 on:
   release:
-    types: [created]
+    types:
+      - created  # 只在 release 创建时触发
 
 jobs:
-  deploy:
-
+  build:
     runs-on: ubuntu-latest
-
+    
     steps:
-    - uses: actions/checkout@v2
+    - name: Checkout repository
+      uses: actions/checkout@v3
+
     - name: Set up Python
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v4
       with:
-        python-version: '3.x'
+        python-version: '>=3.12'  # 可以根据需要修改 Python 版本
+
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
         pip install setuptools wheel twine
-    - name: Build and publish
-      env:
-        TWINE_USERNAME: ${{ secrets.PYPI_USERNAME }}
-        TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }}
+
+    - name: Build distribution
+      run: |
+        python setup.py sdist bdist_wheel  # 使用 setuptools 打包
+
+    - name: Publish to PyPI
       run: |
-        python setup.py sdist bdist_wheel
-        twine upload dist/*
+        twine upload dist/*  # 上传到 PyPI
+      env:
+        TWINE_USERNAME: __token__  # PyPI 推荐使用 '__token__' 作为用户名
+        TWINE_PASSWORD: ${{ secrets.PYPI_API_TOKEN }}  # 使用 GitHub secret 里的 PyPI token
+
diff --git a/.gitignore b/.gitignore
@@ -68,3 +68,4 @@ tags
 .idea
 
 .DS_Store
+.DS_Store
diff --git a/app/fracturex/fracturex/phasefield/main_solve.py b/app/fracturex/fracturex/phasefield/main_solve.py
@@ -243,7 +243,6 @@ def solve_displacement(self) -> float:
         du = self.solver(A, R, atol=1e-20)
         uh += du[:]
         self.uh = uh
-        print('uh', uh)
 
         self.pfcm.update_disp(uh)
         tmr.send('disp_solver')
@@ -294,7 +293,6 @@ def source_coef(bc, index):
         dlform = LinearForm(self.space)
         dlform.add_integrator(ScalarSourceIntegrator(source=source_coef, q=self.q))
         R = dlform.assembly()
-        print('R', R)
         R -= A @ d[:] 
 
 
@@ -304,9 +302,7 @@ def source_coef(bc, index):
         tmr.send('phase_apply_bc')
 
         dd = self.solver(A, R, atol=1e-20)
-        print('R', R)
         d += dd[:]
-        print('d', d)
 
 
         self.d = d

diff --git 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
@@ -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("-----------")