update

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 genarate_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 genarate_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

app/lafem-ims/lafemims/data_generator/near_field_data_generator.py

@@ -2,9 +2,9 @@
 import os
 import matplotlib.pyplot as plt
 from numpy.typing import NDArray
+import numpy as bm
 from typing import Sequence, Callable
 
-from fealpy.backend import backend_manager as bm
 from fealpy.mesh import TriangleMesh, QuadrangleMesh, UniformMesh2d
 from fealpy.pde.pml_2d import PMLPDEModel2d
 from fealpy.functionspace import LagrangeFESpace
@@ -199,7 +199,7 @@ def save(self, save_path: str, scatterer_index: int):
                 d_name = d_values[j]
                 name = f"{k_name}, d={d_name}"
                 data_dict[name] = self.data_for_dsm(k=k_values[i], d=d_values[j])
-        filename = os.path.join(save_path, f"data_for_dsm_{scatterer_index}.bmz")
+        filename = os.path.join(save_path, f"data_for_dsm_{scatterer_index}.npz")
         bm.savez(filename, **data_dict)
 
     def visualization_of_nearfield_data(self, k: float, d: Sequence[float]):

app/lafem-ims/test/test_near_field_data_generator_1.py

@@ -0,0 +1,51 @@
+
+from fealpy.backend import backend_manager as bm
+from fealpy.ml.sampler import CircleCollocator
+
+from lafemims.data_generator import NearFieldDataFEMGenerator2d
+from functional import levelset_circles, genarate_scatterers_circles
+
+
+# 设置随机种子
+SEED = 2024
+
+# 定义计算域
+domain = [-6, 6, -6, 6]
+
+# 定义入射波函数
+u_inc = 'cos(d_0*k*x + d_1*k*y) + sin(d_0*k*x + d_1*k*y) * 1j'
+
+# 定义波矢量方向和波数
+d = [[-bm.sqrt(0.5), bm.sqrt(0.5)]]
+k = [2 * bm.pi]
+
+#散射体个数
+num_of_scatterers = 40000
+# 生成接收点
+num_of_reciever_points = 50
+reciever_points = CircleCollocator(0, 0, 5).run(num_of_reciever_points)
+
+#选择某个散射体
+idx = 0
+centers, radius = genarate_scatterers_circles(num_of_scatterers, SEED)
+
+# 定义水平集函数
+ls_fn = lambda p: levelset_circles(p, centers[idx:idx+1, ...], radius[idx:idx+1, ...])
+
+# 创建近场数据生成器
+generator = NearFieldDataFEMGenerator2d(
+    domain=domain,
+    mesh='UniformMesh',
+    nx=100,
+    ny=100,
+    p=1,
+    q=3,
+    u_inc=u_inc,
+    levelset=ls_fn,
+    d=d,
+    k=k,
+    reciever_points=reciever_points
+)
+
+# 可视化近场数据
+generator.visualization_of_nearfield_data(k=k[0], d=d[0])

app/lafem-ims/test/test_near_field_data_generator_2.py

@@ -0,0 +1,54 @@
+
+from fealpy.backend import backend_manager as bm
+from fealpy.ml.sampler import CircleCollocator
+
+from lafemims.data_generator import NearFieldDataFEMGenerator2d
+from functional import levelset_fourier, genarate_scatterers_fourier
+
+
+# 设置随机种子
+SEED = 2025
+
+# 定义计算域
+domain = [-6, 6, -6, 6]
+
+# 定义入射波函数
+u_inc = 'cos(d_0*k*x + d_1*k*y) + sin(d_0*k*x + d_1*k*y) * 1j'
+
+# 定义波矢量方向和波数
+d = [[-bm.sqrt(0.5), bm.sqrt(0.5)]]
+k = [2 * bm.pi]
+
+#选择散射体复杂度、散射体中心点以及生成散射体个数
+M = 20
+origin_point = bm.array([0.0, 0.0])
+num_of_scatterers = 40000
+
+#确定外层接收点
+num_of_reciever_points = 50
+reciever_points = CircleCollocator(0, 0, 5).run(num_of_reciever_points)
+
+#选择某个散射体
+idx = 0
+coefficients = genarate_scatterers_fourier(M, num_of_scatterers, SEED)     # shape == [2 * M + 1]
+
+# 定义指示函数
+ind_func = lambda p: levelset_fourier(p, M, coefficients[idx, ...], origin_point)
+
+# 创建近场数据生成器
+generator = NearFieldDataFEMGenerator2d(
+    domain=domain,
+    mesh='UniformMesh',
+    nx=100,
+    ny=100,
+    p=1,
+    q=3,
+    u_inc=u_inc,
+    levelset=ind_func,
+    d=d,
+    k=k,
+    reciever_points=reciever_points
+)
+
+# 可视化近场数据
+generator.visualization_of_nearfield_data(k=k[0], d=d[0])

fealpy/fem/scalar_diffusion_integrator.py

@@ -63,7 +63,7 @@ def assembly(self, space: _FS) -> TensorLike:
         bcs, ws, cm = self.fetch(space)
         coef = process_coef_func(coef, bcs=bcs, mesh=mesh, etype='cell', index=self.index)
         gphi = self.fetch_gphix(space)
-
+        
         return bilinear_integral(gphi, gphi, ws, cm, coef, batched=self.batched)
 
     @assemblymethod('fast')