From 289e1152284a05049a8ccc9b1079cc1ca72e2876 Mon Sep 17 00:00:00 2001 From: chaos <2250213115@qq.com> Date: Mon, 25 Nov 2024 19:05:24 +0800 Subject: [PATCH 1/3] update --- app/lafem-eit/script/eit_data_generator.py | 2 - app/lafem-ims/functional.py | 112 ++++++++++ .../near_field_data_generator.py | 4 +- .../default_data/reciever_points_50.npy | Bin 0 -> 928 bytes .../test/test_near_field_data_generator.py | 74 ------- .../test/test_near_field_data_generator_1.py | 51 +++++ .../test/test_near_field_data_generator_2.py | 54 +++++ .../test/test_parallel_data_generation.py | 0 fealpy/fem/scalar_diffusion_integrator.py | 1 - .../ml/generator/near_field_data_generator.py | 200 +++++++++++------- 10 files changed, 347 insertions(+), 151 deletions(-) create mode 100644 app/lafem-ims/functional.py create mode 100644 app/lafem-ims/lafemims/default_data/reciever_points_50.npy delete mode 100644 app/lafem-ims/test/test_near_field_data_generator.py create mode 100644 app/lafem-ims/test/test_near_field_data_generator_1.py create mode 100644 app/lafem-ims/test/test_near_field_data_generator_2.py delete mode 100644 app/lafem-ims/test/test_parallel_data_generation.py diff --git a/app/lafem-eit/script/eit_data_generator.py b/app/lafem-eit/script/eit_data_generator.py index 596c2ff98..e04cdd2f2 100644 --- a/app/lafem-eit/script/eit_data_generator.py +++ b/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) diff --git a/app/lafem-ims/functional.py b/app/lafem-ims/functional.py new file mode 100644 index 000000000..a68703118 --- /dev/null +++ b/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 diff --git a/app/lafem-ims/lafemims/data_generator/near_field_data_generator.py b/app/lafem-ims/lafemims/data_generator/near_field_data_generator.py index 5f8221f79..19bbb4b13 100644 --- a/app/lafem-ims/lafemims/data_generator/near_field_data_generator.py +++ b/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]): diff --git a/app/lafem-ims/lafemims/default_data/reciever_points_50.npy b/app/lafem-ims/lafemims/default_data/reciever_points_50.npy new file mode 100644 index 0000000000000000000000000000000000000000..8e1b75c0a39dd560e56f7f694dbd9d10e7629164 GIT binary patch literal 928 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I%I2098xnmP)#3giMV1~3qDfY1yKi`p+tS#Vw0;l?IT(`CF*?0tl_f>h&$ z9rm>Z+Gs!hY=3UX)VgdE~8U%F*b{MX(>LU7}THbIBWC%;E(xH3Du{9$)cOHaVz z+lEO;UpKKk$fO-_S**e9Aej-WbvTjRLCYy=pG7^V!xJ|n8@@|?4k~7Hr|;&mI*81j zSO4{qfCE$E_Pf8>867HmIbUs)5OP@Z_&3Yy9Y5^lw{o4yU=?-{igc+{{{7b8UTu1} z|59Oxz#TUhoe{rcU%%?z^rv5i9ikxo{V=`)ME!o4`ehLF_QTAB+IOHFYM%zgJqLuL z?$L(0_W;bjvJn3q5QF;X1;oDxzC-;B4UYqtpy6>I5?%+Qpy34#&jYaVgogKl3()X} z#>asgH(>>hA9wb^^2gFgs_kag#_#Y@kjh_R0sPT7TJ~Vz6 zAn|*E5jFk~s6o>QG(HaKLDSD8sC(t1=}QC>9|vAT)8Bijf6>zC0cd NDArray: - """ - 计算水平集函数值。 - - 参数: - - p: 坐标点数组，形状为 (N, 2) - - centers: 圆心坐标数组，形状为 (NCir, 2) - - radius: 半径数组，形状为 (NCir,) - - 返回: - - 水平集函数值，形状为 (N,) - """ - struct = p.shape[:-1] # 获取输入点的结构形状 - p = p.reshape(-1, p.shape[-1]) # 将输入点展平为 (N, 2) - - # 计算每个点到所有圆心的距离 - dis = np.linalg.norm(p[:, None, :] - centers[None, :, :], axis=-1) # 形状为 (N, NCir) - - # 计算每个点到最近圆的距离 - ret = np.min(dis - radius[None, :], axis=-1) # 形状为 (N,) - - return ret.reshape(struct) # 恢复输入点的原始结构形状 - -# 生成接收点 -reciever_points = CircleCollocator(0, 0, 5).run(50) - -# 定义圆的中心和半径 -cirs = np.array([ - [0.4, -0.6, 0.2], - [0.6, -0.5, 0.1], - [0.3, 0.2, 0.3] -], dtype=np.float64) - -centers = cirs[:, 0:2] # 圆心坐标 -radius = cirs[:, 2] # 圆的半径 - -# 定义水平集函数 -ls_fn = lambda p: levelset(p, centers, radius) - -# 创建近场数据生成器 -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]) diff --git a/app/lafem-ims/test/test_near_field_data_generator_1.py b/app/lafem-ims/test/test_near_field_data_generator_1.py new file mode 100644 index 000000000..5437993e3 --- /dev/null +++ b/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]) diff --git a/app/lafem-ims/test/test_near_field_data_generator_2.py b/app/lafem-ims/test/test_near_field_data_generator_2.py new file mode 100644 index 000000000..fe7aeaa97 --- /dev/null +++ b/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]) diff --git a/app/lafem-ims/test/test_parallel_data_generation.py b/app/lafem-ims/test/test_parallel_data_generation.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/fealpy/fem/scalar_diffusion_integrator.py b/fealpy/fem/scalar_diffusion_integrator.py index 4371c3129..f4665c930 100644 --- a/fealpy/fem/scalar_diffusion_integrator.py +++ b/fealpy/fem/scalar_diffusion_integrator.py @@ -63,7 +63,6 @@ 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') diff --git a/fealpy/ml/generator/near_field_data_generator.py b/fealpy/ml/generator/near_field_data_generator.py index 44af20148..92d4f00db 100644 --- a/fealpy/ml/generator/near_field_data_generator.py +++ b/fealpy/ml/generator/near_field_data_generator.py @@ -1,34 +1,40 @@ -import os -from fealpy.backend import backend_manager as bm +import os +import matplotlib.pyplot as plt from numpy.typing import NDArray import numpy as np -from fealpy.solver import spsolve from typing import Sequence, Callable -import matplotlib.pyplot as plt +# from scipy.sparse.linalg import spsolve +# from fealpy.backend import backend_manager as np from fealpy.mesh import TriangleMesh, QuadrangleMesh, UniformMesh2d -from fealpy.functionspace import LagrangeFESpace -from fealpy.fem import ScalarDiffusionIntegrator, ScalarMassIntegrator, ScalarSourceIntegrator, ScalarConvectionIntegrator, DirichletBC -from fealpy.fem import BilinearForm, LinearForm from fealpy.pde.pml_2d import PMLPDEModel2d - -bm.set_backend('numpy') +from fealpy.functionspace import LagrangeFESpace +from fealpy.fem import ( + ScalarDiffusionIntegrator, + ScalarMassIntegrator, + ScalarSourceIntegrator, + ScalarConvectionIntegrator, + BilinearForm, + LinearForm, + DirichletBC +) +from fealpy.solver import spsolve, cg class NearFieldDataFEMGenerator2d: def __init__(self, - domain:Sequence[float], - mesh:str, - nx:int, - ny:int, - p:int, - q:int, - u_inc:str, - levelset:Callable[[NDArray], NDArray], - d:Sequence[float], - k:Sequence[float], - reciever_points:NDArray): + domain: Sequence[float], + mesh: str, + nx: int, + ny: int, + p: int, + q: int, + u_inc: str, + levelset: Callable[[NDArray], NDArray], + d: Sequence[float], + k: Sequence[float], + reciever_points: NDArray): self.domain = domain self.nx = nx @@ -38,52 +44,65 @@ def __init__(self, self.u_inc = u_inc self.levelset = levelset + # 验证并创建网格 if mesh not in ['InterfaceMesh', 'QuadrangleMesh', 'UniformMesh']: raise ValueError("Invalid value for 'mesh'. Choose from 'InterfaceMesh', 'QuadrangleMesh' or 'UniformMesh'.") + + if mesh == 'InterfaceMesh': + self.mesh = TriangleMesh.interfacemesh_generator(box=self.domain, nx=self.nx, ny=self.ny, phi=self.levelset) + self.meshtype = 'InterfaceMesh' + elif mesh == 'QuadrangleMesh': + self.mesh = QuadrangleMesh.from_box(box=self.domain, nx=self.nx, ny=self.ny) + self.meshtype = 'QuadrangleMesh' else: - if mesh == 'InterfaceMesh': - self.mesh = TriangleMesh.interfacemesh_generator(box=self.domain, nx=self.nx, ny=self.ny, phi=self.levelset) - self.meshtype = 'InterfaceMesh' - elif mesh == 'QuadrangleMesh': - self.mesh = QuadrangleMesh.from_box(box=self.domain, nx=self.nx, ny=self.ny) - self.meshtype = 'QuadrangleMesh' - else: - EXTC_1 = self.nx - EXTC_2 = self.ny - HC_1 = 1/EXTC_1 * (self.domain[1] - self.domain[0]) - HC_2 = 1/EXTC_2 * (self.domain[3] - self.domain[2]) - self.mesh = UniformMesh2d((0, EXTC_1, 0, EXTC_2), (HC_1, HC_2), origin=(self.domain[0], self.domain[2])) - self.meshtype = 'UniformMesh' - - self.d = d + EXTC_1 = self.nx + EXTC_2 = self.ny + HC_1 = 1 / EXTC_1 * (self.domain[1] - self.domain[0]) + HC_2 = 1 / EXTC_2 * (self.domain[3] - self.domain[2]) + self.mesh = UniformMesh2d((0, EXTC_1, 0, EXTC_2), (HC_1, HC_2), origin=(self.domain[0], self.domain[2])) + self.meshtype = 'UniformMesh' + + self.d = d self.k = k self.reciever_points = reciever_points qf = self.mesh.quadrature_formula(self.q, 'cell') - self.bc, _= qf.get_quadrature_points_and_weights() - - def get_nearfield_data(self, k:float, d:Sequence[float]): - - k_index = (self.k).index(k) - d_index = (self.d).index(d) - pde = PMLPDEModel2d(levelset=self.levelset, - domain=self.domain, - u_inc=self.u_inc, - A=1, - k=self.k[k_index], - d=self.d[d_index], - refractive_index=[1, 1+1/self.k[k_index]**2], - absortion_constant=1.79, - lx=1.0, - ly=1.0 - ) + self.bc, _ = qf.get_quadrature_points_and_weights() + + def get_nearfield_data(self, k: float, d: Sequence[float]): + """ + 获取近场数据。 + + 参数: + - k: 波数 + - d: 波矢量方向 + + 返回: + - uh: 近场数据 + """ + k_index = self.k.index(k) + d_index = self.d.index(d) + pde = PMLPDEModel2d( + levelset=self.levelset, + domain=self.domain, + u_inc=self.u_inc, + A=1, + k=self.k[k_index], + d=self.d[d_index], + refractive_index=[1, 1 + 1 / self.k[k_index]**2], + absortion_constant=1.79, + lx=1.0, + ly=1.0 + ) space = LagrangeFESpace(self.mesh, p=self.p) + # 定义积分器 D = ScalarDiffusionIntegrator(pde.diffusion_coefficient, q=self.q) C = ScalarConvectionIntegrator(pde.convection_coefficient, q=self.q) M = ScalarMassIntegrator(pde.reaction_coefficient, q=self.q) f = ScalarSourceIntegrator(pde.source, q=self.q) + # 组装双线性形式和线性形式 b = BilinearForm(space) b.add_integrator([D, C, M]) @@ -92,14 +111,30 @@ def get_nearfield_data(self, k:float, d:Sequence[float]): A = b.assembly() F = l.assembly() + bc = DirichletBC(space, pde.dirichlet) - uh = space.function(dtype=bm.complex128) + uh = space.function(dtype=np.complex128) A, F = bc.apply(A, F) uh[:] = spsolve(A, F, solver='scipy') + # uh[:] = cg(A, F) + # uh[:] = spsolve(A.to_scipy(), F) + print(uh.shape) return uh - - def points_location_and_bc(self, p, domain:Sequence[float], nx:int, ny:int): + def points_location_and_bc(self, p: NDArray, domain: Sequence[float], nx: int, ny: int): + """ + 计算接收点的位置和重心坐标。 + + 参数: + - p: 接收点坐标 + - domain: 计算域 + - nx: x方向的网格数 + - ny: y方向的网格数 + + 返回: + - location: 接收点所在单元的索引 + - bc: 重心坐标 + """ x = p[..., 0] y = p[..., 1] cell_length_x = (domain[1] - domain[0]) / nx @@ -110,22 +145,31 @@ def points_location_and_bc(self, p, domain:Sequence[float], nx:int, ny:int): bc_x_ = ((x - domain[0]) / cell_length_x) % 1 bc_y_ = ((y - domain[2]) / cell_length_y) % 1 - bc_x = bm.array([[bc_x_, 1 - bc_x_]], dtype=bm.float64) - bc_y = bm.array([[bc_y_, 1 - bc_y_]], dtype=bm.float64) + bc_x = np.array([[bc_x_, 1 - bc_x_]], dtype=np.float64) + bc_y = np.array([[bc_y_, 1 - bc_y_]], dtype=np.float64) bc = (bc_x, bc_y) return location, bc - def data_for_dsm(self, k:float, d:Sequence[float]): + def data_for_dsm(self, k: float, d: Sequence[float]): + """ + 获取用于DSM的数据。 + + 参数: + - k: 波数 + - d: 波矢量方向 + 返回: + - data: DSM数据 + """ reciever_points = self.reciever_points data_length = reciever_points.shape[0] - data = bm.zeros((data_length, ), dtype=bm.complex128) + data = np.zeros((data_length,), dtype=np.complex128) uh = self.get_nearfield_data(k=k, d=d) - - if self.meshtype =='InterfaceMesh': + + if self.meshtype == 'InterfaceMesh': b = self.mesh.point_to_bc(reciever_points) location = self.mesh.location(reciever_points) - for i in range (data_length): + for i in range(data_length): data[i] = uh(b[i])[location[i]] elif self.meshtype == 'QuadrangleMesh': for i in range(data_length): @@ -142,14 +186,20 @@ def data_for_dsm(self, k:float, d:Sequence[float]): u = uh(b).reshape(-1) data[i] = u[location] return data - - def save(self, save_path:str, scatterer_index:int): + def save(self, save_path: str, scatterer_index: int): + """ + 保存数据。 + + 参数: + - save_path: 保存路径 + - scatterer_index: 散射体索引 + """ k_values = self.k d_values = self.d data_dict = {} - for i in range (len(k_values)): - for j in range (len(d_values)): + for i in range(len(k_values)): + for j in range(len(d_values)): k_name = f'k={k_values[i]}' d_name = d_values[j] name = f"{k_name}, d={d_name}" @@ -157,24 +207,30 @@ def save(self, save_path:str, scatterer_index:int): filename = os.path.join(save_path, f"data_for_dsm_{scatterer_index}.npz") np.savez(filename, **data_dict) - def visualization_of_nearfield_data(self, k:float, d:Sequence[float]): + def visualization_of_nearfield_data(self, k: float, d: Sequence[float]): + """ + 可视化近场数据。 + 参数: + - k: 波数 + - d: 波矢量方向 + """ uh = self.get_nearfield_data(k=k, d=d) value = uh(self.bc) if self.meshtype == 'UniformMesh': - self.mesh.ftype = bm.float64 + self.mesh.ftype = np.float64 self.mesh.add_plot(plt, cellcolor=value[..., 0].real, linewidths=0) self.mesh.add_plot(plt, cellcolor=value[..., 0].imag, linewidths=0) - #TODO + # TODO: 添加更多可视化选项 # fig = plt.figure() # axes = fig.add_subplot(1, 3, 1) # self.mesh.add_plot(axes) # if self.meshtype == 'UniformMesh': - # uh = uh.view(bm.ndarray) + # uh = uh.view(np.ndarray) # axes = fig.add_subplot(1, 3, 2, projection='3d') - # self.mesh.show_function(axes, bm.real(uh)) + # self.mesh.show_function(axes, np.real(uh)) # axes = fig.add_subplot(1, 3, 3, projection='3d') - # self.mesh.show_function(axes, bm.imag(uh)) + # self.mesh.show_function(axes, np.imag(uh)) plt.show() \ No newline at end of file From 4c85bd38c7696655ec7e24e606da3710dcae3e4b Mon Sep 17 00:00:00 2001 From: chaos <2250213115@qq.com> Date: Mon, 25 Nov 2024 19:07:57 +0800 Subject: [PATCH 2/3] - --- fealpy/fem/scalar_diffusion_integrator.py | 1 + 1 file changed, 1 insertion(+) diff --git a/fealpy/fem/scalar_diffusion_integrator.py b/fealpy/fem/scalar_diffusion_integrator.py index f4665c930..3b21c1dad 100644 --- a/fealpy/fem/scalar_diffusion_integrator.py +++ b/fealpy/fem/scalar_diffusion_integrator.py @@ -63,6 +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') From cb98c39c6a6660039513abcc8859a915d9c86897 Mon Sep 17 00:00:00 2001 From: chaos <2250213115@qq.com> Date: Mon, 25 Nov 2024 19:09:53 +0800 Subject: [PATCH 3/3] - --- .../ml/generator/near_field_data_generator.py | 27 ++++++++----------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/fealpy/ml/generator/near_field_data_generator.py b/fealpy/ml/generator/near_field_data_generator.py index 92d4f00db..19bbb4b13 100644 --- a/fealpy/ml/generator/near_field_data_generator.py +++ b/fealpy/ml/generator/near_field_data_generator.py @@ -2,11 +2,9 @@ import os import matplotlib.pyplot as plt from numpy.typing import NDArray -import numpy as np +import numpy as bm from typing import Sequence, Callable -# from scipy.sparse.linalg import spsolve -# from fealpy.backend import backend_manager as np from fealpy.mesh import TriangleMesh, QuadrangleMesh, UniformMesh2d from fealpy.pde.pml_2d import PMLPDEModel2d from fealpy.functionspace import LagrangeFESpace @@ -19,7 +17,7 @@ LinearForm, DirichletBC ) -from fealpy.solver import spsolve, cg +from fealpy.solver import spsolve class NearFieldDataFEMGenerator2d: @@ -113,12 +111,9 @@ def get_nearfield_data(self, k: float, d: Sequence[float]): F = l.assembly() bc = DirichletBC(space, pde.dirichlet) - uh = space.function(dtype=np.complex128) + uh = space.function(dtype=bm.complex128) A, F = bc.apply(A, F) uh[:] = spsolve(A, F, solver='scipy') - # uh[:] = cg(A, F) - # uh[:] = spsolve(A.to_scipy(), F) - print(uh.shape) return uh def points_location_and_bc(self, p: NDArray, domain: Sequence[float], nx: int, ny: int): @@ -145,8 +140,8 @@ def points_location_and_bc(self, p: NDArray, domain: Sequence[float], nx: int, n bc_x_ = ((x - domain[0]) / cell_length_x) % 1 bc_y_ = ((y - domain[2]) / cell_length_y) % 1 - bc_x = np.array([[bc_x_, 1 - bc_x_]], dtype=np.float64) - bc_y = np.array([[bc_y_, 1 - bc_y_]], dtype=np.float64) + bc_x = bm.array([[bc_x_, 1 - bc_x_]], dtype=bm.float64) + bc_y = bm.array([[bc_y_, 1 - bc_y_]], dtype=bm.float64) bc = (bc_x, bc_y) return location, bc @@ -163,7 +158,7 @@ def data_for_dsm(self, k: float, d: Sequence[float]): """ reciever_points = self.reciever_points data_length = reciever_points.shape[0] - data = np.zeros((data_length,), dtype=np.complex128) + data = bm.zeros((data_length,), dtype=bm.complex128) uh = self.get_nearfield_data(k=k, d=d) if self.meshtype == 'InterfaceMesh': @@ -205,7 +200,7 @@ def save(self, save_path: str, scatterer_index: int): 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}.npz") - np.savez(filename, **data_dict) + bm.savez(filename, **data_dict) def visualization_of_nearfield_data(self, k: float, d: Sequence[float]): """ @@ -218,7 +213,7 @@ def visualization_of_nearfield_data(self, k: float, d: Sequence[float]): uh = self.get_nearfield_data(k=k, d=d) value = uh(self.bc) if self.meshtype == 'UniformMesh': - self.mesh.ftype = np.float64 + self.mesh.ftype = bm.float64 self.mesh.add_plot(plt, cellcolor=value[..., 0].real, linewidths=0) self.mesh.add_plot(plt, cellcolor=value[..., 0].imag, linewidths=0) @@ -227,10 +222,10 @@ def visualization_of_nearfield_data(self, k: float, d: Sequence[float]): # axes = fig.add_subplot(1, 3, 1) # self.mesh.add_plot(axes) # if self.meshtype == 'UniformMesh': - # uh = uh.view(np.ndarray) + # uh = uh.view(bm.ndarray) # axes = fig.add_subplot(1, 3, 2, projection='3d') - # self.mesh.show_function(axes, np.real(uh)) + # self.mesh.show_function(axes, bm.real(uh)) # axes = fig.add_subplot(1, 3, 3, projection='3d') - # self.mesh.show_function(axes, np.imag(uh)) + # self.mesh.show_function(axes, bm.imag(uh)) plt.show() \ No newline at end of file