Merge remote-tracking branch 'upstream/master'

app/FuelRodSim/HeatEquationData.py

@@ -113,7 +113,7 @@ def domain(self):
     def duration(self):
         return [0, 1]
 
-    def source(self,p,t):
+    def source(self, p, t):
         return 0
 
     def dirichlet(self, p, t):

app/lafem-eit/lafemeit/model/data_feature.py

@@ -13,11 +13,6 @@
 class DataPreprocessor(nn.Module):
     """Data Feature Boundary Solver based on FEM."""
     def __init__(self, solver: LaplaceFEMSolver) -> None:
-        """_summary_
-
-        Args:
-            solver (LaplaceFEMSolver): _description_
-        """
         super().__init__()
         self.solver = solver
         self.vuh = None
@@ -33,7 +28,6 @@ def forward(self, input: Tensor) -> Tensor:
         gd, gn = input[:, 0, :], input[:, 1, :] # [B*CH, NN_bd]
         vuh = solver.solve_from_potential(gd) # [B*CH, gdof]
         vn = solver.normal_derivative(vuh) # [B*CH, bddof]
-        # gnvn = solver.solve_from_gnf(None, gn-vn, f_only=True)[:, solver.bd_dof_flag] # [B*CH, bddof]
         gnvn = gn - vn # [B*CH, bddof]
         self.vuh = vuh
 
@@ -44,12 +38,6 @@ class DataFeature(nn.Module):
     """Data Feature Solver based on FEM."""
     def __init__(self, solver: LaplaceFEMSolver,
                  bc_filter: Optional[TensorFunc]=None) -> None:
-        """_summary_
-
-        Args:
-            solver (LaplaceFEMSolver): _description_
-            bc_filter (Optional[TensorFunc], optional): _description_. Defaults to None.
-        """
         super().__init__()
         self.solver = solver
         self.bc_filter = bc_filter

app/lafem-eit/lafemeit/model/unet100.py

@@ -136,7 +136,6 @@ def build_eit_model(
     ):
     BdDofs = ext * 4
     mesh = TriangleMesh.from_box([-1, 1, -1, 1], ext, ext, device=device)
-    print(mesh.boundary_face_index().shape)
 
     if fractype == 'nograd':
         frac = Fractional(BdDofs, device=device)

app/lafem-eit/script/eit_data_generator.py

@@ -19,7 +19,6 @@
 from fealpy.cem import EITDataGenerator
 
 
-logger.setLevel('INFO')
 parser = argparse.ArgumentParser()
 parser.add_argument("config", help="path of the .yaml file")
 
@@ -73,9 +72,10 @@ def neumann(points: Tensor, *args):
     return bm.sin(bm.tensordot(freq, theta, axes=0))
 
 
-def main(sigma_iterable: Sequence[int], seed=0, index=0):
+def main(sigma_iterable: Sequence[int], seed: int = 0, index: int = 0):
     np.random.seed(seed)
     umesh = UniformMesh2d([0, EXT, 0, EXT], [H, H], [-1., -1.], itype=bm.int32, ftype=DTYPE)
+    pixel = umesh.entity('node')
 
     for sigma_idx in tqdm(sigma_iterable,
                           desc=f"task{index}",
@@ -96,7 +96,7 @@ def main(sigma_iterable: Sequence[int], seed=0, index=0):
         interface_mesh = TriangleMesh.interfacemesh_generator(umesh, phi=ls_fn)
         generator = EITDataGenerator(mesh=interface_mesh, p=P, q=Q)
         gn = generator.set_boundary(neumann, batch_size=len(FREQ))
-        label = generator.set_levelset(SIGMA, ls_fn)
+        label = generator.set_levelset(SIGMA, ls_fn, pixel)
         gd = generator.run()
 
         np.save(

app/lafem-eit/script/laplace_beltrami_eigen.py

@@ -34,7 +34,7 @@
 P_ = config['fem']['order']
 
 
-def laplace_eigen_fem(mesh: IntervalMesh, p: int=1) -> Tuple[NDArray, NDArray, NDArray, NDArray]:
+def laplace_eigen_fem(mesh: IntervalMesh, p: int = 1) -> Tuple[NDArray, NDArray, NDArray, NDArray]:
     """_summary_
 
     Args:
@@ -58,6 +58,7 @@ def laplace_eigen_fem(mesh: IntervalMesh, p: int=1) -> Tuple[NDArray, NDArray, N
 EXTx, EXTy = config['mesh']["ext"]
 Lx, Ly = config['mesh']['length']
 Origin = config['mesh']['origin']
+refine = config['mesh']['refine']
 box = [Origin[0], Origin[0] + Lx, Origin[1], Origin[1] + Ly]
 
 print("Generating Boundary Laplace Beltrami operator...")
@@ -75,12 +76,19 @@ def laplace_eigen_fem(mesh: IntervalMesh, p: int=1) -> Tuple[NDArray, NDArray, N
     mesh = IntervalMesh.from_mesh_boundary(tri_mesh)
     del tri_mesh
 
+    bform_2 = BilinearForm(LagrangeFESpace(mesh))
+    bform_2.add_integrator(ScalarMassIntegrator(q=Q_))
+    M0 = bform_2.assembly().to_dense()
+    del bform_2
+
     NN = mesh.number_of_nodes()
+    mesh.uniform_refine(refine)
 
-    w, v, _, M = laplace_eigen_fem(mesh, p=P_)
+    w, v, _, _ = laplace_eigen_fem(mesh, p=P_)
     w = w[1:NN+1]
-    vinv = (v.T @ M)[1:NN+1, :NN]
+    vinv = v.T[1:NN+1, :NN] @ M0
     v = v[:NN, 1:NN+1]
+    print(vinv @ v)
 
     np.savez(config['file'], w=w, v=v, vinv=vinv)
     print("Saved.")

example/experiment/fem/maxwell_eq_n2fem_dirichlet.py

@@ -67,7 +67,7 @@
 maxit = args.maxit
 maxit = 3
 dim = args.dim
-dim = 3
+#dim = 3
 backend = args.backend
 bm.set_backend(backend)
 if dim == 2:

example/opt/Levy_example.py

@@ -26,8 +26,8 @@
     levy_steps_tensor = bm.array(levy_steps, device=device, dtype=bm.float32)
     bround_steps_tensor = bm.array(bround_steps, device=device, dtype=bm.float32)
 
-    L[n * i : (i + 1) * n, :, :] = bm.cumsum(levy_steps_tensor, axis=0)
-    B[n * i : (i + 1) * n, :, :] = bm.cumsum(bround_steps_tensor, axis=0)
+    L[1 : 1 + n + 1, :, :] = bm.cumsum(levy_steps_tensor, axis=0)
+    B[1 : 1 + n + 1, :, :] = bm.cumsum(bround_steps_tensor, axis=0)
 
 end_time = time.perf_counter()
 running_time = end_time - start_time
@@ -38,15 +38,38 @@
 L_cpu = L.cpu().numpy()
 B_cpu = B.cpu().numpy()
 
-for particle in range(0, 1):
-    L_x_vals = L_cpu[:, 0, particle]  
-    L_y_vals = L_cpu[:, 1, particle]  
-    B_x_vals = B_cpu[:, 0, particle]  
-    B_y_vals = B_cpu[:, 1, particle]  
+
+plt.scatter(L_cpu[0, 0, 0], L_cpu[0, 0, 0], marker=MarkerStyle('*'), color='red', s=500) 
+plt.xticks([])
+plt.yticks([])
+
+for i in range(1, 16):
+    if i < 6:
+        step = 50
+    elif i < 11:
+        step = 150
+    else:
+        step = 300
+    L_x_vals = L_cpu[ : step * i, 0, 0]
+    L_y_vals = L_cpu[ : step * i, 1, 0]
+    B_x_vals = B_cpu[ : step * i, 0, 0]
+    B_y_vals = B_cpu[ : step * i, 1, 0] 
     plt.figure()
     plt.plot(L_x_vals, L_y_vals, '-', color = 'blue')
     plt.plot(B_x_vals, B_y_vals, '-', color = 'black')
     plt.scatter(L_x_vals[0], L_y_vals[0], marker=MarkerStyle('*'), color='red', s=500) 
+    plt.xticks([])
+    plt.yticks([])
+
+particle = 0
+L_x_vals = L_cpu[ : , 0, particle]
+L_y_vals = L_cpu[ : , 1, particle]
+B_x_vals = B_cpu[ : , 0, particle]
+B_y_vals = B_cpu[ : , 1, particle] 
+plt.figure()
+plt.plot(L_x_vals, L_y_vals, '-', color = 'blue')
+plt.plot(B_x_vals, B_y_vals, '-', color = 'black')
+plt.scatter(L_x_vals[0], L_y_vals[0], marker=MarkerStyle('*'), color='red', s=500) 
 plt.xticks([])
 plt.yticks([])
 plt.show()

example/opt/MultiAGV_example.py

@@ -9,7 +9,8 @@
 import sys
 import time 
 from fealpy.backend import backend_manager as bm
-from fealpy.iopt.A_star import AStar, GridMap, Graph
+from fealpy.opt.A_star import AStar, GridMap, Graph
+# from fealpy.iopt.A_star import AStar, GridMap, Graph
 bm.set_backend('pytorch')
 
 

example/opt/ant_example.py

@@ -1,7 +1,7 @@
 import time
 from fealpy.backend import backend_manager as bm
 import matplotlib.pyplot as plt
-from fealpy.iopt.ANT_TSP import calD, Ant_TSP
+from fealpy.opt.ANT_TSP import calD, Ant_TSP
 # bm.set_backend('pytorch')
 
 # 导入数据(34个城市)

example/opt/pso_example.py

@@ -1,4 +1,4 @@
-from fealpy.iopt.particle_swarm_opt_alg import  PSOProblem, PSO, QPSO
+from fealpy.opt.particle_swarm_opt_alg import PSOProblem, PSO, QPSO
 from fealpy.backend import backend_manager as bm
 import time
 

fealpy/cem/generator/eit_data_generator.py

@@ -23,8 +23,8 @@ def __init__(self, mesh: Mesh, p: int=1, q: Optional[int]=None) -> None:
         Args:
             mesh (Mesh): _description_
             p (int, optional): Order of the Lagrange finite element space. Defaults to 1.
-            q (int | None, optional): Order of the quadrature, use `q = p + 2` if None.\
-            Defaults to None.
+            q (int | None, optional): Order of the quadrature, use `q = p + 2` if None.
+                Defaults to None.
         """
         q = p + 2 if q is None else q
 
@@ -53,16 +53,19 @@ def __init__(self, mesh: Mesh, p: int=1, q: Optional[int]=None) -> None:
         self.cdata_indices = bm.stack([cdata_row, cdata_col], axis=0)
 
     def set_levelset(self, sigma_vals: Tuple[float, float],
-                     levelset: Callable[[Tensor], Tensor]) -> Tensor:
+                     levelset: Callable[[Tensor], Tensor],
+                     pixel: Tensor) -> Tensor:
         """Set inclusion distribution.
 
         Args:
             sigma_vals (Tuple[float, float]): Sigma value of inclusion and background.
             levelset (Callable): level-set function indicating inclusion and background.
+            pixel (Tensor): Sampler points on the interior on which the inclusion is indicated.
+                This may be different from the mesh nodes or cell barycenters.
 
         Returns:
             Tensor: Label boolean Tensor with True for inclusion nodes and False\
-                for background nodes, shaped (nodes, ).
+                for background nodes, shaped (pixels, ).
         """
         def _coef_func(p: Tensor):
             inclusion = levelset(p) < 0.
@@ -78,17 +81,16 @@ def _coef_func(p: Tensor):
         self._di.coef = _coef_func
         self._di.clear() # clear the cached result as the coef has changed
         bform.add_integrator(self._di)
-        self._A = bform.assembly()
+        self._A = bform.assembly(format='coo')
 
         cdata_indices = self.cdata_indices
         cdataT_indices = bm.flip(cdata_indices, axis=0)
         A_n_indices = bm.concat([self._A.indices(), cdata_indices, cdataT_indices], axis=1)
         A_n_values = bm.concat([self._A.values(), self.cdata, self.cdata], axis=-1)
-        self.A_n = COOTensor(A_n_indices, A_n_values, spshape=(self.gdof+1, self.gdof+1))
-        self.A_n = self.A_n.tocsr()
+        A_n = COOTensor(A_n_indices, A_n_values, spshape=(self.gdof+1, self.gdof+1))
+        self.A_n = A_n.tocsr()
 
-        node = space.mesh.entity('node')
-        return levelset(node) < 0.
+        return levelset(pixel) < 0.
 
     def set_boundary(self, gn_source: Union[Callable[[Tensor], Tensor], Tensor],
                      batch_size: int=0, *, zero_integral=False) -> Tensor:

fealpy/fem/bilinear_form.py

@@ -74,7 +74,7 @@ def _scalar_assembly(self, retain_ints: bool, batch_size: int):
         return M
 
     @overload
-    def assembly(self, *, retain_ints: bool=False) -> COOTensor: ...
+    def assembly(self, *, retain_ints: bool=False) -> CSRTensor: ...
     @overload
     def assembly(self, *, format: Literal['coo'], retain_ints: bool=False) -> COOTensor: ...
     @overload

fealpy/fem/dirichlet_bc.py

@@ -230,7 +230,7 @@ def apply_vector(self, vector: TensorLike, matrix: SparseTensor,
         else:
             if uh is None:
                 uh = bm.zeros_like(f)
-            uh = self.space.boundary_interpolate(gd, uh, self.threshold)
+            uh, _ = self.space.boundary_interpolate(gd, uh, self.threshold)
 
         bd_idx = self.boundary_dof_index
         f = f - A.matmul(uh[:])

fealpy/fem/mthlaplace_integrator.py

@@ -76,8 +76,6 @@ def assembly_without_numerical_integration(self, space: _FS):
         q = space.p+3 if self.q is None else self.q
         cmcoeff = space.coeff
         bgm = MLaplaceBernsteinIntegrator(m=m, q=q).assembly(space.bspace)
-        print(bgm.dtype)
-        print(cmcoeff.dtype)
         M = bm.einsum('cil,clm,cpm->cip', cmcoeff, bgm, cmcoeff)
         return M
 

fealpy/fem/vector_mass_integrator.py

@@ -68,14 +68,15 @@ def assembly_cell_matrix_for_vector_basis_vspace(self, space, index=_S, cellmeas
         coef = self.coef
         mesh = space.mesh
         GD = mesh.geo_dimension()
+        self.device = mesh.device
 
         if cellmeasure is None:
             cellmeasure = mesh.entity_measure('cell', index=index)
 
         NC = len(cellmeasure)
         ldof = space.number_of_local_dofs() 
         if out is None:
-            D = bm.zeros((NC, ldof, ldof), dtype=space.ftype)
+            D = bm.zeros((NC, ldof, ldof), device=self.device,dtype=space.ftype)
         else:
             D = out
 

fealpy/functionspace/bernstein_fe_space.py

@@ -20,6 +20,7 @@ class BernsteinFESpace(FunctionSpace, Generic[_MT]):
     def __init__(self, mesh: _MT, p: int=1, ctype='C'):
         self.mesh = mesh
         self.p = p
+        self.device = mesh.device
 
         assert ctype in {'C', 'D'}
         self.ctype = ctype # 空间连续性类型
@@ -85,19 +86,18 @@ def basis(self, bcs: TensorLike, index: Index=_S, p = None):
         ldof = multiIndex.shape[0]
 
         B = bcs
-        B = bm.ones((NQ, p+1, TD+1), dtype=self.ftype)
+        B = bm.ones((NQ, p+1, TD+1), device=self.device, dtype=self.ftype)
         # B[:, 1:] = bcs[:, None, :]
         B = bm.set_at(B,(slice(None),slice(1,None)),bcs[:,None,:])
         B = bm.cumprod(B, axis=1)
-        P = bm.arange(p+1)
+        P = bm.arange(p+1, device=self.device)
         P = bm.set_at(P,(0),1)
         P = bm.cumprod(P,axis=0).reshape(1, -1, 1)
         B = B/P
 
         # B : (NQ, p+1, TD+1) 
         # B[:, multiIndex, bm.arange(TD+1).reshape(1, -1)]: (NQ, ldof, TD+1)
-        phi = P[0, -1, 0]*bm.prod(B[:, multiIndex, 
-            bm.arange(TD+1).reshape(1, -1)], axis=-1)
+        phi = P[0, -1, 0]*bm.prod(B[:, multiIndex, bm.arange(TD+1).reshape(1, -1)], axis=-1)
         return phi[None, :]
 
     @barycentric
@@ -132,26 +132,26 @@ def grad_basis(self, bcs: TensorLike, index: Index=_S, variable='u',p=None):
         ldof = multiIndex.shape[0]
 
         B = bcs
-        B = bm.ones((NQ, p+1, TD+1), dtype=self.ftype)
+        B = bm.ones((NQ, p+1, TD+1), device=self.device, dtype=self.ftype)
         # B[:, 1:] = bcs[:, None, :]
         B = bm.set_at(B,(slice(None),slice(1,None)),bcs[:,None,:])
         B = bm.cumprod(B, axis=1)
 
-        P = bm.arange(p+1)
+        P = bm.arange(p+1,device=self.device)
         # P[0] = 1
         P = bm.set_at(P,(0),1)
         P = bm.cumprod(P,axis=0).reshape(1, -1, 1)
         B = B/P
 
-        F = bm.zeros(B.shape, dtype=self.ftype)
+        F = bm.zeros(B.shape, device=self.device, dtype=self.ftype)
         # F[:, 1:] = B[:, :-1]
         F = bm.set_at(F,(slice(None),slice(1,None)),B[:,:-1])
         shape = bcs.shape[:-1]+(ldof, TD+1)
-        R = bm.zeros(shape, dtype=self.ftype)
+        R = bm.zeros(shape,device=self.device,dtype=self.ftype)
         for i in range(TD+1):
             idx = list(range(TD+1))
             idx.remove(i)
-            idx = bm.array(idx, dtype=self.itype)
+            idx = bm.array(idx,device=self.device, dtype=self.itype)
             # R[..., i] = bm.prod(B[..., multiIndex[:, idx], idx.reshape(1, -1)],axis=-1)*F[..., multiIndex[:, i], [i]]
             R = bm.set_at(R,(...,i),bm.prod(B[..., multiIndex[:, idx], idx.reshape(1, -1)],axis=-1)*F[..., multiIndex[:, i], [i]])
         Dlambda = self.mesh.grad_lambda()
@@ -167,7 +167,7 @@ def hess_basis(self, bcs: TensorLike, index: Index=_S, variable='u'):
         g2phi = self.grad_m_basis(bcs, 2, index=index)
         TD = self.mesh.top_dimension()
         shape = g2phi.shape[:-1] + (TD, TD)
-        hval  = bm.zeros(shape, dtype=self.ftype)
+        hval  = bm.zeros(shape, device=self.device, dtype=self.ftype)
         hval = bm.set_at(hval,(...,0,0),g2phi[..., 0])
         hval = bm.set_at(hval,(...,0,1),g2phi[..., 1])
         hval = bm.set_at(hval,(...,1,0),g2phi[..., 1])
@@ -222,7 +222,7 @@ def grad_m_basis(self, bcs: TensorLike, m: int, index = _S):
         midxp_1 = bm.multi_index_matrix(m, GD) # m   次多重指标
 
         N, N1 = len(symidx), midxp_1.shape[0]
-        B = bm.zeros((N1, NQ, ldof), dtype=self.ftype)
+        B = bm.zeros((N1, NQ, ldof), device=self.device, dtype=self.ftype)
         symLambdaBeta = bm.zeros((N1, NC, N), dtype=self.ftype)
         for beta, Bi, symi in zip(midxp_1, B, symLambdaBeta):
             midxp_0 -= beta[None, :]
@@ -336,7 +336,7 @@ def boundary_interpolate(self, gD: Union[Callable, int, float],
         isbdface = self.mesh.boundary_face_flag()
         f2d  = self.face_to_dof()[isbdface]
 
-        isDDof = bm.zeros(gdof, dtype=bm.bool)
+        isDDof = bm.zeros(gdof, device=self.device, dtype=bm.bool)
         # isDDof[f2d] = True
         isDDof = bm.set_at(isDDof,(f2d),True)
         if callable(gD):