Merge remote-tracking branch 'upstream/develop' into develop

.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
+

.gitignore

@@ -68,3 +68,4 @@ tags
 .idea
 
 .DS_Store
+.DS_Store

app/fracturex/fracturex/cases/phase_field/model3d.py

@@ -68,8 +68,8 @@ def is_dirchlet_boundary(self, p):
 
 
 parser.add_argument('--mesh_type',
-        default='tri', type=str,
-        help='网格类型, 默认为 tri.')
+        default='tet', type=str,
+        help='网格类型, 默认为 tet.')
 
 parser.add_argument('--enable_adaptive',
         default=False, type=bool,
@@ -166,8 +166,8 @@ def is_dirchlet_boundary(self, p):
 tmr.send(None)
 end = time.time()
 
-force = ms.Rforce
-disp = ms.force_value
+force = ms.get_residual_force()
+disp = model.is_z_force()
 
 ftname = 'force_'+args.mesh_type + '_p' + str(p) + '_' + 'model3d_disp.pt'
 
@@ -177,6 +177,7 @@ def is_dirchlet_boundary(self, p):
 with open(tname, 'w') as file:
     file.write(f'\n time: {end-start},\n degree:{p},\n, backend:{backend},\n, model_type:{model_type},\n, enable_adaptive:{enable_adaptive},\n, marking_strategy:{marking_strategy},\n, refine_method:{refine_method},\n, n:{n},\n, maxit:{maxit},\n, vtkname:{vtkname}\n')
 fig, axs = plt.subplots()
+disp = model.is_z_force()
 plt.plot(disp, force, label='Force')
 plt.xlabel('Displacement Increment')
 plt.ylabel('Residual Force')

app/fracturex/fracturex/cases/phase_field/square_domian_with_fracture.py

@@ -19,7 +19,7 @@ def __init__(self):
         E = 210
         nu = 0.3
         Gc = 2.7e-3
-        l0 = 0.02
+        l0 = 0.015
         self.params = {'E': E, 'nu': nu, 'Gc': Gc, 'l0': l0}
 
     def is_y_force(self):
@@ -204,16 +204,25 @@ def adaptive_mesh(self, mesh, d0=0.49, d1=1.01, h=0.005):
 tmr.send(None)
 end = time.time()
 
-force = ms.Rforce
-disp = ms.force_value
+force = ms.get_residual_force()
+
 
 ftname = 'force_'+args.mesh_type + '_p' + str(p) + '_' + 'model1_disp.pt'
 
 torch.save(force, ftname)
 #np.savetxt('force'+tname, bm.to_numpy(force))
+
 tname = 'params_'+args.mesh_type + '_p' + str(p) + '_' + 'model1_disp.txt'
 with open(tname, 'w') as file:
     file.write(f'\n time: {end-start},\n degree:{p},\n, backend:{backend},\n, model_type:{model_type},\n, enable_adaptive:{enable_adaptive},\n, marking_strategy:{marking_strategy},\n, refine_method:{refine_method},\n, n:{n},\n, maxit:{maxit},\n, vtkname:{vtkname}\n')
+
+if force_type == 'y':
+    disp = model.is_y_force()
+elif force_type == 'x':
+    disp = model.is_x_force()
+else:
+    raise ValueError('Invalid force type.')
+
 fig, axs = plt.subplots()
 plt.plot(disp, force, label='Force')
 plt.xlabel('Displacement Increment')

app/fracturex/fracturex/phasefield/main_solve.py

@@ -75,7 +75,7 @@ def __init__(self, mesh, material_params: Dict,
         self.tmr = timer()
         next(self.tmr)
 
-    def initialization_settings(self, p: int = 1, q: int = None, ):
+    def initialize_settings(self, p: int = 1, q: int = None, ):
         """
         Initialize the settings for the problem.
         """
@@ -121,11 +121,11 @@ def solve(self, method: str = 'lfem', p: int = 1, q: int = None, maxit: int = 50
             VTK output file name, by default None.
         """
         self._method = method
-        self.initialization_settings(p=p, q=q)
+        self.initialize_settings(p=p, q=q)
         self._initialize_force_boundary()
         self._Rforce = bm.zeros_like(self._force_value)
 
-#        for i in range(1):
+        #for i in range(2):
         for i in range(len(self._force_value)-1):
             print('i', i)
             self._currt_force_value = self._force_value[i+1]
@@ -258,31 +258,40 @@ def solve_phase_field(self) -> float:
             The norm of the residual.
         """
         Gc, l0, d = self.Gc, self.l0, self.d
+
+        @barycentric
+        def diff_coef(bc, index):
+            gg_hd, c_d = self.CSDFunc.grad_grad_density_function(d(bc))
+            return Gc * l0 * 2 / c_d
+
+        @barycentric
+        def mass_coef1(bc, index):
+            gg_hd, c_d = self.CSDFunc.grad_grad_density_function(d(bc))
+            return gg_hd * Gc / (l0 * c_d)
+
 
         @barycentric
-        def coef(bc, index):
-            gg_gd = self.EDFunc.grad_grad_degradation_function(d)
+        def mass_coef2(bc, index):
+            gg_gd = self.EDFunc.grad_grad_degradation_function(d(bc))
             return gg_gd * self.pfcm.maximum_historical_field(bc)
 
         @barycentric
         def source_coef(bc, index):
-            gg_gd = self.EDFunc.grad_degradation_function_constant_coef()
-            return -1 * gg_gd * self.pfcm.maximum_historical_field(bc)
+            gc_gd = self.EDFunc.grad_degradation_function_constant_coef()
+            return -1 * gc_gd * self.pfcm.maximum_historical_field(bc)
 
         tmr = self.tmr
         tmr.send('phase_start')
 
-        gg_hd, c_d = self.CSDFunc.grad_grad_density_function(d)
-
         dbform = BilinearForm(self.space)
-        dbform.add_integrator(ScalarDiffusionIntegrator(coef=Gc * l0 * 2 / c_d, q=self.q))
-        dbform.add_integrator(ScalarMassIntegrator(coef=gg_hd * Gc / (l0 * c_d), q=self.q))
-        dbform.add_integrator(ScalarMassIntegrator(coef=source_coef, q=self.q))
+        dbform.add_integrator(ScalarDiffusionIntegrator(coef=diff_coef, q=self.q))
+        dbform.add_integrator(ScalarMassIntegrator(coef=mass_coef1, q=self.q))
+        dbform.add_integrator(ScalarMassIntegrator(coef=mass_coef2, q=self.q))
         A = dbform.assembly()
         tmr.send('phase_matrix_assemble')
 
         dlform = LinearForm(self.space)
-        dlform.add_integrator(ScalarSourceIntegrator(source=coef, q=self.q))
+        dlform.add_integrator(ScalarSourceIntegrator(source=source_coef, q=self.q))
         R = dlform.assembly()
         R -= A @ d[:] 
 
@@ -409,11 +418,6 @@ def mass_kernel_func2(u):
         def mass_grad_kernel_func2(u):
             return self.EDFunc.grad_grad_degradation_function(u)
 
-        @barycentric
-        def source_coef(bc, index):
-            gg_gd = self.EDFunc.grad_degradation_function_constant_coef()
-            return -1 * gg_gd * self.pfcm.maximum_historical_field(bc)
-
         diffusion_coef.kernel_func = diffusion_kernel_func
         mass_coef1.kernel_func = mass_kernel_func1
         mass_coef2.kernel_func = mass_kernel_func2
@@ -667,7 +671,7 @@ def set_energy_degradation(self, degradation_type='quadratic', EDfunc=None, **kw
             Additional parameters passed to the energy degradation function.
         """
         if EDfunc is not None:
-            self.EDFunc = EDfunc(degradation_type=degradation_type, **kwargs)
+            self.EDFunc = EDfunc(degradation_type='user_defined', **kwargs)
         else:
             self.EDFunc = EDFunc(degradation_type=degradation_type)
 
@@ -685,7 +689,7 @@ def set_crack_surface_density(self, density_type='AT2', CSDfunc=None, **kwargs):
             Additional parameters passed to the crack surface density function.
         """
         if CSDfunc is not None:
-            self.CSDFunc = CSDfunc(density_type=density_type, **kwargs)
+            self.CSDFunc = CSDfunc(density_type=='user_defined', **kwargs)
         else:
             self.CSDFunc = CSDFunc(density_type=density_type)
 

app/fracturex/fracturex/phasefield/phase_fracture_material.py

@@ -1,4 +1,5 @@
 import numpy as np
+import torch
 from typing import Optional
 
 from fealpy.typing import TensorLike
@@ -286,7 +287,18 @@ def strain_pm_eig_decomposition(self, s: TensorLike):
         
         @param[in] s strain，（NC, NQ, GD, GD）
         """
-        w, v = bm.linalg.eigh(s) # w 特征值, v 特征向量
+        '''
+        if bm.device_type(s) == 'cuda':
+            torch.cuda.empty_cache()
+            try:
+                w, v = bm.linalg.eigh(s)  # w 特征值, v 特征向量
+            except torch.cuda.OutOfMemoryError as e:
+                print("CUDA out of memory. Attempting to free cache.")
+                torch.cuda.empty_cache()
+        else:
+            w, v = bm.linalg.eigh(s) # w 特征值, v 特征向量
+        '''
+        w, v = bm.linalg.eigh(s)
         p, m = self.macaulay_operation(w)
 
         sp = bm.zeros_like(s)
@@ -322,6 +334,16 @@ def heaviside(self, x):
         val[x < -1e-13] = 0
         return val
 
+    def linear_strain_value(self, bc):
+        """
+        Compute the linear strain tensor.
+        """
+        bc = bm.array([[1/3, 1/3, 1/3]], dtype=bm.float64)
+        uh = self.uh
+        guh = uh.grad_value(bc)
+        strain = 0.5 * (guh + bm.swapaxes(guh, -2, -1))
+        return strain
+
     @ barycentric
     def maximum_historical_field(self, bc):
 

app/lafem-eit/script/eit_data_generator.py

@@ -104,9 +104,7 @@ def main(sigma_iterable: Sequence[int], seed: int = 0, index: int = 0):
         ctrs_ = np.random.rand(NUM_CIR, 2) * 1.6 - 0.8 # (NCir, GD)
         b = np.min(0.9-np.abs(ctrs_), axis=-1) # (NCir, )
         rads_ = np.random.rand(NUM_CIR) * (b-0.1) + 0.1 # (NCir, )
-        ctrs = bm.from_numpy(ctrs_)
         ctrs = bm.astype(ctrs, bm.float64)
-        rads = bm.from_numpy(rads_)
         rads = bm.astype(rads_, bm.float64)
 
         ls_fn = lambda p: levelset(p, ctrs, rads)

app/lafem-ims/functional.py

@@ -0,0 +1,112 @@
+
+from numpy.typing import NDArray
+from fealpy.backend import backend_manager as bm
+from fealpy.typing import TensorLike
+
+#多圆形并集-型 散射体
+def levelset_circles(p: NDArray, centers: NDArray, radius: NDArray) -> NDArray:
+    """
+    参数:
+    - p: 坐标点数组，形状为 (N, 2)
+    - centers: 圆心坐标数组，形状为 (NCir, 2)
+    - radius: 半径数组，形状为 (NCir,)
+    
+    返回:
+    - 水平集函数值，形状为 (N,)
+    """
+    struct = p.shape[:-1]  # 获取输入点的结构形状
+    p = p.reshape(-1, p.shape[-1])  # 将输入点展平为 (N, 2)
+
+    # 计算每个点到所有圆心的距离
+    dis = bm.linalg.norm(p[:, None, :] - centers[None, :, :], axis=-1)  # 形状为 (N, NCir)
+
+    # 计算每个点到最近圆的距离
+    ret = bm.min(dis - radius[None, :], axis=-1)  # 形状为 (N,)
+
+    return ret.reshape(struct)  # 恢复输入点的原始结构形状
+
+#半径函数为傅里叶系数-型 散射体
+def levelset_fourier(points: NDArray, complexity: int, coefficient: NDArray, origin: NDArray) -> NDArray:
+    """
+    参数:
+    - points: NDArray, 形状为 (N, 2)，表示 N 个点的坐标。
+    - complexity: int, 水平集函数的复杂度。
+    - coefficient: NDArray, 形状为 (2 * M + 1,)，表示傅里叶系数。
+    - origin: NDArray, 形状为 (2,)，表示原点坐标。
+
+    返回:
+    - flag: NDArray, 形状为 (N,)，表示每个点是否在水平集函数内。
+    """
+    points = points - origin
+    angles_rad = bm.zeros_like(points[:, 0], dtype=bm.float64)  # 创建一个和点集大小相同的数组，用于存储角度弧度
+
+    # 处理分母为零的情况
+    zero_indices = bm.where(points[:, 0] == 0)
+    angles_rad[zero_indices] = bm.pi / 2 if bm.any(points[zero_indices, 1] > 0) else 3 * bm.pi / 2
+
+    # 处理分母不为零的情况
+    non_zero_indices = bm.where(points[:, 0] != 0)
+    slopes = points[non_zero_indices, 1] / points[non_zero_indices, 0]  # 计算斜率
+    angles_rad[non_zero_indices] = bm.arctan(slopes)  # 计算角度弧度
+
+    # 将负值转换为正值
+    negative_angle_indices = bm.where(angles_rad < 0)
+    angles_rad[negative_angle_indices] += bm.pi
+
+    # 调整角度弧度，确保在0到2*pi之间
+    angles_rad = angles_rad % (2 * bm.pi)
+
+    # 处理负斜率的情况
+    negative_slope_indices = bm.where((points[:, 0] >= 0) & (points[:, 1] < 0))
+    angles_rad[negative_slope_indices] += bm.pi
+
+    r_t = coefficient[0]
+
+    for i in range(complexity):
+        r_t += coefficient[i + 1] * bm.cos((i + 1) * angles_rad) + coefficient[i + complexity + 1] * bm.sin((i + 1) * angles_rad)
+
+    distances = r_t - bm.linalg.norm(points, axis=1)
+    flag = distances >= 0
+
+    return flag
+
+#################################################################################################################################
+
+def generate_scatterers_circles(num_of_scatterers: int, seed: int) ->TensorLike:
+    """
+    参数:
+    - num_of_scatterers: int, 散射体的数量。
+    
+    返回:
+    - ctrs: TensorLike, 形状为 (num_of_scatterers, 2)，表示每个散射体的中心坐标。
+    - rads: TensorLike, 形状为 (num_of_scatterers,)，表示每个散射体的半径。
+    """
+    bm.random.seed(seed)
+
+    ctrs = bm.random.rand(num_of_scatterers, 2) * 1.6 - 0.8 # (NCir, GD)
+    b = bm.min(0.9-bm.abs(ctrs), axis=-1) # (NCir, )
+    rads = bm.random.rand(num_of_scatterers) * (b-0.1) + 0.1 # (NCir, )
+    ctrs = bm.astype(ctrs, bm.float64)
+    rads = bm.astype(rads, bm.float64)
+
+    return ctrs, rads
+def generate_scatterers_fourier(complexity: int, num_of_scatterers: int, seed: int) ->TensorLike:
+    """
+    参数:
+    - complexity: int, 水平集函数的复杂度。
+    - num_of_scatterers: int, 散射体的数量。
+    
+    返回:
+    - c: TensorLike, 形状为 (num_of_scatterers, 2*complexity+1)，表示每个散射体的参数。
+    """
+    bm.random.seed(seed)
+
+    c = bm.zeros((num_of_scatterers, 2*complexity+1), dtype=bm.float64)
+    radius = bm.random.uniform(0, 0.1, (num_of_scatterers, complexity))
+    theta = bm.random.uniform(0, 2*bm.pi, (num_of_scatterers, complexity))
+
+    c[:, 0:1] = bm.random.uniform(1, 1.2, (num_of_scatterers, 1))
+    c[:, 1:complexity+1] = radius * bm.cos(theta)
+    c[:, complexity+1: ] = radius * bm.sin(theta)
+
+    return c