From 3b926f1881715aa4c1cf75372a28f95faeea8e68 Mon Sep 17 00:00:00 2001 From: wpx <1143615697@qq.com> Date: Fri, 18 Oct 2024 09:20:27 +0800 Subject: [PATCH] update --- .../fem/doublelaplace_fem_dirichlet.py | 20 +- .../fem/navier_stokes_cylinder_2d.py | 202 +++++------------- .../fem/navier_stokes_cylinder_2d_new.py | 139 ------------ fealpy/fem/mthlaplace_integrator.py | 4 +- fealpy/functionspace/bernstein_fe_space.py | 2 +- .../functionspace/cm_conforming_fe_space.py | 19 +- fealpy/functionspace/functional.py | 12 +- fealpy/mesh/triangle_mesh.py | 11 +- fealpy/pde/biharmonic_triharmonic_2d.py | 8 +- 9 files changed, 106 insertions(+), 311 deletions(-) delete mode 100644 example/experiment/fem/navier_stokes_cylinder_2d_new.py diff --git a/example/experiment/fem/doublelaplace_fem_dirichlet.py b/example/experiment/fem/doublelaplace_fem_dirichlet.py index 958d20231..871daf224 100644 --- a/example/experiment/fem/doublelaplace_fem_dirichlet.py +++ b/example/experiment/fem/doublelaplace_fem_dirichlet.py @@ -34,11 +34,11 @@ help='初始网格剖分段数.') parser.add_argument('--maxit', - default=4, type=int, + default=2, type=int, help='默认网格加密求解的次数, 默认加密求解 4 次') parser.add_argument('--backend', - default='numpy', type=str, + default='pytorch', type=str, help='默认后端为numpy') parser.add_argument('--meshtype', @@ -49,7 +49,7 @@ bm.set_backend(args.backend) -decive = "cuda" +device = "cuda" p = args.degree n = args.n meshtype = args.meshtype @@ -62,19 +62,23 @@ u = (sp.sin(2*sp.pi*y)*sp.sin(2*sp.pi*x))**2 pde = DoubleLaplacePDE(u) ulist = get_flist(u)[:3] -mesh = TriangleMesh.from_box([0,1,0,1], n, n) -NDof = bm.zeros(maxit, dtype=bm.float64) +mesh = TriangleMesh.from_box([0,1,0,1], n, n, device=device) + +ikwargs = bm.context(mesh.cell) +fkwargs = bm.context(mesh.node) + +NDof = bm.zeros(maxit, **ikwargs) errorType = ['$|| u - u_h||_{\\Omega,0}$', '$||\\nabla u - \\nabla u_h||_{\\Omega,0}$', '$||\\nabla^2 u - \\nabla^2 u_h||_{\\Omega,0}$'] -errorMatrix = bm.zeros((3, maxit), dtype=bm.float64) +errorMatrix = bm.zeros((3, maxit), **fkwargs) tmr.send('网格和pde生成时间') for i in range(maxit): node = mesh.entity('node') - isCornerNode = bm.zeros(len(node),dtype=bm.bool) - for n in bm.array([[0,0],[1,0],[0,1],[1,1]], dtype=bm.float64): + isCornerNode = bm.zeros(len(node),dtype=bm.bool, device=device) + for n in bm.array([[0,0],[1,0],[0,1],[1,1]], **fkwargs): isCornerNode = isCornerNode | (bm.linalg.norm(node-n[None, :], axis=1)<1e-10) diff --git a/example/experiment/fem/navier_stokes_cylinder_2d.py b/example/experiment/fem/navier_stokes_cylinder_2d.py index e80336369..9a676ec42 100644 --- a/example/experiment/fem/navier_stokes_cylinder_2d.py +++ b/example/experiment/fem/navier_stokes_cylinder_2d.py @@ -7,8 +7,8 @@ @bref @ref ''' -from fealpy import logger -logger.setLevel('ERROR') +#from fealpy import logger +#logger.setLevel('ERROR') from fealpy.backend import backend_manager as bm from fealpy.mesh import TriangleMesh @@ -16,14 +16,14 @@ from fealpy.functionspace import TensorFunctionSpace from fealpy.fem import ( BilinearForm, ScalarDiffusionIntegrator, - ScalarMassIntegrator, PressWorkIntegrator0, PressWorkIntegrator , - PressWorkIntegrator1, ScalarConvectionIntegrator) + ScalarMassIntegrator, PressWorkIntegrator , + ScalarConvectionIntegrator) from fealpy.fem import LinearForm, ScalarSourceIntegrator from fealpy.fem import DirichletBC from fealpy.sparse import COOTensor from fealpy.fem import VectorSourceIntegrator -from fealpy.fem import BlockForm +from fealpy.fem import LinearBlockForm, BlockForm from fealpy.solver import spsolve from fealpy.pde.navier_stokes_equation_2d import FlowPastCylinder @@ -45,58 +45,30 @@ rho = pde.rho mu = pde.mu -mesh = pde.mesh1(0.05) - +mesh = pde.mesh(0.05,'fealpy') timeline = UniformTimeLine(0, T, nt) dt = timeline.dt pspace = LagrangeFESpace(mesh, p=pdegree) -uspace = LagrangeFESpace(mesh, p=udegree) -tensor_uspace = TensorFunctionSpace(uspace, (2, -1)) - -u0x = uspace.function() -u0y = uspace.function() -u0 = bm.stack((u0x[:], u0y[:]), axis=1) -u1x = uspace.function() -u1y = uspace.function() +space = LagrangeFESpace(mesh, p=udegree) +uspace = TensorFunctionSpace(space, (2, -1)) + +u0 = uspace.function() +u1 = uspace.function() p1 = pspace.function() ugdof = uspace.number_of_global_dofs() pgdof = pspace.number_of_global_dofs() -gdof = pgdof+2*ugdof +gdof = ugdof + pgdof fname = output + 'test_'+ str(0).zfill(10) + '.vtu' -mesh.nodedata['velocity'] = u0 +mesh.nodedata['velocity'] = u1.reshape(2,-1).T mesh.nodedata['pressure'] = p1 mesh.to_vtk(fname=fname) -APX_bform = BilinearForm((pspace, uspace)) -APX_bform.add_integrator(PressWorkIntegrator0(coef=-1, q=q)) - -APY_bform = BilinearForm((pspace, uspace)) -APY_bform.add_integrator(PressWorkIntegrator1(coef=-1, q=q)) -''' -M_bform = BilinearForm(uspace) -M_bform.add_integrator(ScalarMassIntegrator(rho/dt, q=q)) -M = M_bform.assembly() - -S_bform = BilinearForm(uspace) -S_bform.add_integrator(ScalarDiffusionIntegrator(mu, q=q)) -S = S_bform.assembly() - -P_bform = BilinearForm((pspace, tensor_uspace)) +##BilinearForm +P_bform = BilinearForm((pspace, uspace)) P_bform.add_integrator(PressWorkIntegrator(-1, q=q)) -P = P_bform.assembly() - - -from scipy.sparse import bmat -A = bmat([[APX.to_scipy()],[APY.to_scipy()]], format='coo') - -print(bm.sum(bm.abs(A.toarray()-P.to_scipy().toarray()))) - -exit() -''' - A_bform = BilinearForm(uspace) A_bform.add_integrator(ScalarMassIntegrator(rho/dt, q=q)) @@ -104,128 +76,70 @@ ConvectionIntegrator = ScalarConvectionIntegrator(q=q) A_bform.add_integrator(ConvectionIntegrator) -APX_bform = BilinearForm((pspace, uspace)) -APX_bform.add_integrator(PressWorkIntegrator0(coef=-1, q=q)) -APX = APX_bform.assembly() - -APY_bform = BilinearForm((pspace, uspace)) -APY_bform.add_integrator(PressWorkIntegrator1(coef=-1, q=q)) -APY = APY_bform.assembly() +##LinearForm +ulform = LinearForm(uspace) +SourceIntegrator = ScalarSourceIntegrator(q = q) +ulform.add_integrator(SourceIntegrator) +SourceIntegrator.source = u0 +plform = LinearForm(pspace) +b = LinearBlockForm([ulform, plform]) +bb = b.assembly() +print(bb) +exit() #边界处理 +## 边界处理太繁琐 +## threshold一定要传tuple或者list吗 +u_isbddof = uspace.is_boundary_dof() +u_out_isbd = uspace.is_boundary_dof(threshold=(pde.is_outflow_boundary,)) +u_isbddof[u_out_isbd] = False +p_isbddof = pspace.is_boundary_dof(threshold=pde.is_outflow_boundary) +isBdDof = bm.concatenate([u_isbddof, p_isbddof]) +## b +#xu,u_in_isbd = uspace.boundary_interpolate(pde.u_inflow_dirichlet, threshold=(pde.is_inflow_boundary,)) +#xp = bm.zeros(pgdof) +#axx = bm.concatenate((xu,xp)) + xx = bm.zeros(gdof, dtype=mesh.ftype) -u_isbddof_u0 = uspace.is_boundary_dof() -u_isbddof_in = uspace.is_boundary_dof(threshold = pde.is_inflow_boundary) -u_isbddof_out = uspace.is_boundary_dof(threshold = pde.is_outflow_boundary) -u_isbddof_circle = uspace.is_boundary_dof(threshold = pde.is_circle_boundary) - -u_isbddof_u0[u_isbddof_in] = False -u_isbddof_u0[u_isbddof_out] = False -xx[0:ugdof][u_isbddof_u0] = 0 -xx[ugdof:2*ugdof][u_isbddof_u0] = 0 - -u_isbddof = u_isbddof_u0 -u_isbddof[u_isbddof_in] = True -#ipoint = uspace.interpolation_points()[u_isbddof_in] -ipoint = uspace.interpolation_points() +u_isbddof_in = space.is_boundary_dof(threshold = pde.is_inflow_boundary) +ipoint = space.interpolation_points() uinflow = pde.u_inflow_dirichlet(ipoint) - -#xx[0:ugdof][u_isbddof_in] = uinflow[:,0] -#xx[ugdof:2*ugdof][u_isbddof_in] = uinflow[:,1] -#print(u_bd_in.shape) -#print(xx.shape) -#xx[u_bd_in] = uinflow.reshape(-1) - -p_isBdDof_p0 = pspace.is_boundary_dof(threshold = pde.is_outflow_boundary) -bd = bm.concatenate((u_isbddof_in, u_isbddof_in, p_isBdDof_p0)) +bd = bm.concatenate((u_isbddof_in, u_isbddof_in, p_isbddof)) value_bd = bm.concatenate((uinflow[:,0],uinflow[:,1], bm.zeros(pgdof))) xx[bd] = value_bd[bd] -isBdDof = bm.concatenate([u_isbddof, u_isbddof, p_isBdDof_p0], axis=0) - for i in range(10): t1 = timeline.next_time_level() print("time=", t1) - - @barycentric - def concoef(bcs, index): - a1 = u0x(bcs, index) - a2 = u0y(bcs, index) - result = bm.concatenate((a1[...,bm.newaxis],a2[..., bm.newaxis]), axis=2) - return result - - ConvectionIntegrator.coef = concoef + ConvectionIntegrator.coef = u0 ConvectionIntegrator.clear() - A = BlockForm([[A_bform, None, APX_bform], - [None, A_bform, APY_bform], - [APX_bform.T, APY_bform.T, None]]) + A = BlockForm([[A_bform, P_bform], + [P_bform.T, None]]) A = A.assembly() - ''' - if backend == 'numpy': - from scipy.sparse import coo_array, bmat - def coo(A): - data = A._values - indices = A._indices - return coo_array((data, indices)) - - A = bmat([[coo(M+S+C), None,coo(-APX)], - [None, coo(M+S+C), coo(-APY)], - [coo(-APX).T, coo(-APY).T, None]], format='coo') - A = COOTensor(bm.stack([A.row,A.col],axis=0), A.data, spshape=A.shape) - A = A.coalesce() - if backend == 'pytorch': - indices = bm.tensor([[],[]]) - data = bm.tensor([]) - zeros_0 = COOTensor(indices, data, (ugdof,ugdof)) - zeros_1 = COOTensor(indices, data, (pgdof,pgdof)) - A0 = bm.concatenate([M+S+C, zeros_0, -APX], axis=1) - A1 = bm.concatenate([zeros_0, M+S+C, -APY], axis=1) - A2 = bm.concatenate([-APX.T, -APY.T ,zeros_1], axis=1) - A = bm.concatenate((A0,A1,A2),axis=0) - ''' - lform = LinearForm(uspace) - lform.add_integrator(ScalarSourceIntegrator(u0x)) - b0 = lform.assembly() - lform = LinearForm(uspace) - lform.add_integrator(ScalarSourceIntegrator(u0y)) - b1 = lform.assembly() + SourceIntegrator.source = u0 + SourceIntegrator.clear() + ##LinearBlockForm + b0 = ulform.assembly() b2 = bm.zeros(pgdof) - b = bm.concatenate([b0,b1,b2]) - - b -= A@xx - b[isBdDof] = xx[isBdDof] + b = bm.concatenate([b0, b2]) - A = DirichletBC(uspace, xx, isDDof=isBdDof).apply_matrix(A, check=False) - #A,b = DirichletBC(uspace, xx, isDDof=isBdDof).apply(A, b, check=None) - ''' - import scipy.sparse as sp - values = A.values() - indices = A.indices() - A = sp.coo_matrix((values, (indices[0], indices[1])), shape=A.shape) - A = A.tocsr() - x = sp.linalg.spsolve(A,b) - ''' + A,b = DirichletBC((uspace,pspace), xx, threshold=isBdDof).apply(A, b) x = spsolve(A, b, 'mumps') - x = bm.array(x) - ''' - x = cg(A, b, maxiter=10000) - ''' - u1x[:] = x[:ugdof] - u1y[:] = x[ugdof:2*ugdof] - p1[:] = x[2*ugdof:] - - u0x[:] = u1x[:] - u0y[:] = u1y[:] - u0 = bm.stack((u0x[:], u0y[:]), axis=1) + u1[:] = x[:ugdof] + p1[:] = x[ugdof:] + fname = output + 'test_'+ str(i+1).zfill(10) + '.vtu' - mesh.nodedata['velocity'] = u0 + ##mesh.nodedata + mesh.nodedata['velocity'] = u1.reshape(2,-1).T mesh.nodedata['pressure'] = p1 mesh.to_vtk(fname=fname) - - timeline.advance() + + u0[:] = u1[:] + timeline.advance() +print(bm.sum(bm.abs(u1))) diff --git a/example/experiment/fem/navier_stokes_cylinder_2d_new.py b/example/experiment/fem/navier_stokes_cylinder_2d_new.py deleted file mode 100644 index 1f806c399..000000000 --- a/example/experiment/fem/navier_stokes_cylinder_2d_new.py +++ /dev/null @@ -1,139 +0,0 @@ -#!/usr/bin/python3 -'''! - @Author: wpx - @File Name: navier_stokes_cylinder_2d.py - @Mail: wpx15673207315@gmail.com - @Created Time: Mon 12 Aug 2024 04:52:25 PM CST - @bref - @ref -''' -#from fealpy import logger -#logger.setLevel('ERROR') - -from fealpy.backend import backend_manager as bm -from fealpy.mesh import TriangleMesh -from fealpy.functionspace import LagrangeFESpace -from fealpy.functionspace import TensorFunctionSpace -from fealpy.fem import ( - BilinearForm, ScalarDiffusionIntegrator, - ScalarMassIntegrator, PressWorkIntegrator0, PressWorkIntegrator , - PressWorkIntegrator1, ScalarConvectionIntegrator) - -from fealpy.fem import LinearForm, ScalarSourceIntegrator -from fealpy.fem import DirichletBC -from fealpy.sparse import COOTensor -from fealpy.fem import VectorSourceIntegrator -from fealpy.fem import BlockForm -from fealpy.solver import spsolve - -from fealpy.pde.navier_stokes_equation_2d import FlowPastCylinder -from fealpy.decorator import barycentric, cartesian -from fealpy.old.timeintegratoralg import UniformTimeLine -from fealpy.fem import DirichletBC - -#backend = 'pytorch' -backend = 'numpy' -bm.set_backend(backend) - -output = './' -udegree = 2 -pdegree = 1 -q = 4 -T = 5 -nt = 5000 -pde = FlowPastCylinder() -rho = pde.rho -mu = pde.mu - -mesh = pde.mesh(0.05,'fealpy') -timeline = UniformTimeLine(0, T, nt) -dt = timeline.dt - -pspace = LagrangeFESpace(mesh, p=pdegree) -space = LagrangeFESpace(mesh, p=udegree) -uspace = TensorFunctionSpace(space, (2, -1)) - -u0 = uspace.function() -u1 = uspace.function() -p1 = pspace.function() - -ugdof = uspace.number_of_global_dofs() -pgdof = pspace.number_of_global_dofs() -gdof = ugdof + pgdof - -fname = output + 'test_'+ str(0).zfill(10) + '.vtu' -mesh.nodedata['velocity'] = u1.reshape(2,-1).T -mesh.nodedata['pressure'] = p1 -mesh.to_vtk(fname=fname) - -##BilinearForm -P_bform = BilinearForm((pspace, uspace)) -P_bform.add_integrator(PressWorkIntegrator(-1, q=q)) - -A_bform = BilinearForm(uspace) -A_bform.add_integrator(ScalarMassIntegrator(rho/dt, q=q)) -A_bform.add_integrator(ScalarDiffusionIntegrator(mu, q=q)) -ConvectionIntegrator = ScalarConvectionIntegrator(q=q) -A_bform.add_integrator(ConvectionIntegrator) - -##LinearForm -lform = LinearForm(uspace) -SourceIntegrator = ScalarSourceIntegrator(q = q) -lform.add_integrator(SourceIntegrator) - - -#边界处理 -## 边界处理太繁琐 -## threshold一定要传tuple或者list吗 -u_isbddof = uspace.is_boundary_dof() -u_out_isbd = uspace.is_boundary_dof(threshold=(pde.is_outflow_boundary,)) -u_isbddof[u_out_isbd] = False -p_isbddof = pspace.is_boundary_dof(threshold=pde.is_outflow_boundary) -isBdDof = bm.concatenate([u_isbddof, p_isbddof]) -## b -#xu,u_in_isbd = uspace.boundary_interpolate(pde.u_inflow_dirichlet, threshold=(pde.is_inflow_boundary,)) -#xp = bm.zeros(pgdof) -#axx = bm.concatenate((xu,xp)) - -xx = bm.zeros(gdof, dtype=mesh.ftype) -u_isbddof_in = space.is_boundary_dof(threshold = pde.is_inflow_boundary) -ipoint = space.interpolation_points() -uinflow = pde.u_inflow_dirichlet(ipoint) -bd = bm.concatenate((u_isbddof_in, u_isbddof_in, p_isbddof)) -value_bd = bm.concatenate((uinflow[:,0],uinflow[:,1], bm.zeros(pgdof))) -xx[bd] = value_bd[bd] - - -for i in range(10): - t1 = timeline.next_time_level() - print("time=", t1) - - ConvectionIntegrator.coef = u0 - ConvectionIntegrator.clear() - - A = BlockForm([[A_bform, P_bform], - [P_bform.T, None]]) - A = A.assembly() - - SourceIntegrator.source = u0 - SourceIntegrator.clear() - ##LinearBlockForm - b0 = lform.assembly() - b2 = bm.zeros(pgdof) - b = bm.concatenate([b0, b2]) - - A,b = DirichletBC((uspace,pspace), xx, threshold=isBdDof).apply(A, b) - x = spsolve(A, b, 'mumps') - - u1[:] = x[:ugdof] - p1[:] = x[ugdof:] - - - fname = output + 'test_'+ str(i+1).zfill(10) + '.vtu' - ##mesh.nodedata - mesh.nodedata['velocity'] = u1.reshape(2,-1).T - mesh.nodedata['pressure'] = p1 - mesh.to_vtk(fname=fname) - - u0[:] = u1[:] - timeline.advance() diff --git a/fealpy/fem/mthlaplace_integrator.py b/fealpy/fem/mthlaplace_integrator.py index a324ee9ee..60a92273f 100644 --- a/fealpy/fem/mthlaplace_integrator.py +++ b/fealpy/fem/mthlaplace_integrator.py @@ -57,9 +57,11 @@ def assembly(self, space: _FS) -> TensorLike: m = self.m coef = self.coef mesh = getattr(space, 'mesh', None) + device = bm.get_device(mesh.cell) GD = mesh.geo_dimension() idx = bm.array(mesh.multi_index_matrix(m, GD-1)) - num = factorial(m)/bm.prod(bm.array(factorial(idx)), axis=1) + idx_cpu = bm.to_numpy(idx) + num = factorial(m)/bm.prod(bm.array(factorial(idx_cpu), device=device), axis=1) bcs, ws, gmphi, cm, index = self.fetch(space) coef = process_coef_func(coef, bcs=bcs, mesh=mesh, etype='cell', index=index) diff --git a/fealpy/functionspace/bernstein_fe_space.py b/fealpy/functionspace/bernstein_fe_space.py index 3da0e6fd5..812be5381 100644 --- a/fealpy/functionspace/bernstein_fe_space.py +++ b/fealpy/functionspace/bernstein_fe_space.py @@ -223,7 +223,7 @@ def grad_m_basis(self, bcs: TensorLike, m: int, index = _S): N, N1 = len(symidx), midxp_1.shape[0] B = bm.zeros((N1, NQ, ldof), device=self.device, dtype=self.ftype) - symLambdaBeta = bm.zeros((N1, NC, N), dtype=self.ftype) + symLambdaBeta = bm.zeros((N1, NC, N), dtype=self.ftype, device=self.device) for beta, Bi, symi in zip(midxp_1, B, symLambdaBeta): midxp_0 -= beta[None, :] idx = bm.where(bm.all(midxp_0>-1, axis=1))[0] diff --git a/fealpy/functionspace/cm_conforming_fe_space.py b/fealpy/functionspace/cm_conforming_fe_space.py index c9d3591fb..1bc27e740 100644 --- a/fealpy/functionspace/cm_conforming_fe_space.py +++ b/fealpy/functionspace/cm_conforming_fe_space.py @@ -159,6 +159,7 @@ def coefficient_matrix(self): p = self.p m = self.m mesh = self.mesh + device = bm.get_device(mesh.cell) isCornerNode = self.isCornerNode NC = mesh.number_of_cells() @@ -201,9 +202,13 @@ def coefficient_matrix(self): j = midx2num(beta[None, :]) j = dof2num[j] dbeta = beta[flag] - coeff[:, i, j] = sign*comb(bm.sum(beta), - r)*bm.prod(bm.array(comb(dalpha, dbeta)), - axis=0)*factorial(r) + beta_cpu = bm.to_numpy(bm.sum(beta)) + dalpha_cpu = bm.to_numpy(dalpha) + dbeta_cpu = bm.to_numpy(dbeta) + coeff[:, i, j] = sign*comb(beta_cpu,r)*bm.prod(bm.array( + comb(dalpha_cpu, dbeta_cpu),device=device),axis=0)*factorial(r) + #coeff[:, i, j] = sign*comb(bm.sum(beta),r)*bm.prod(bm.array(comb(dalpha, dbeta),**ikwargs), + # axis=0)*factorial(r) # 局部自由度 边 @@ -221,11 +226,13 @@ def coefficient_matrix(self): dalpha = alpha[de] betas = mesh.multi_index_matrix(int(dalpha), 2) for beta in betas: - val = bm.prod(N[:, de]**beta,axis=1)/bm.prod(bm.array(factorial(beta)), axis=0) + beta_cpu = bm.to_numpy(beta) + val = bm.prod(N[:, de]**beta,axis=1)/bm.prod(bm.array(factorial(beta_cpu), device=device), axis=0) beta[e] = beta[e] +alpha[e] + beta_sum_cpu = bm.to_numpy(bm.sum(beta)) j = midx2num(beta[None, :]) j = dof2num[j] - coeff[:, i, j] = comb(bm.sum(beta), + coeff[:, i, j] = comb(beta_sum_cpu, int(dalpha))*(factorial(int(dalpha)))**2*val[:, None] #for i in range(NC): @@ -247,7 +254,7 @@ def coefficient_matrix(self): Ndelat[:, 2, 1] = -t0 ndof = self.number_of_internal_dofs('node') coeff1 = bm.zeros((NC, 3*ndof, 3*ndof), dtype=self.ftype, device=self.device) - symidx = [symmetry_index(2, r) for r in range(1, 2*m+1)] + symidx = [symmetry_index(2, r, device=device) for r in range(1, 2*m+1)] # 边界自由度 isBdNode = mesh.boundary_node_flag() isBdEdge = mesh.boundary_face_flag() diff --git a/fealpy/functionspace/functional.py b/fealpy/functionspace/functional.py index 680a50504..dff45f2b0 100644 --- a/fealpy/functionspace/functional.py +++ b/fealpy/functionspace/functional.py @@ -123,15 +123,17 @@ def symmetry_span_array(arr, alpha): break return ret -def symmetry_index(d, r, dtype=None): +def symmetry_index(d, r, dtype=None, device=None): dtype = dtype if dtype is not None else bm.int32 """ @brief 将 d 维 r 阶张量拉长以后，其对称部分对应的索引和出现的次数 """ symidx0 = bm.tensor(list(combinations_with_replacement(range(d), r)), - dtype=dtype) - coe = bm.flip(d**bm.arange(r, dtype=dtype)) - symidx = bm.einsum('ij,j->i', symidx0, coe) + dtype=dtype, device=device) + coe = bm.flip(d**bm.arange(r, dtype=dtype, device=device)) + + symidx = bm.einsum('ij,j->i', bm.astype(symidx0, bm.float64), bm.astype(coe, bm.float64)) + symidx = bm.astype(symidx, dtype) midx = bm.multi_index_matrix(r, d-1) #midx0 = bm.zeros_like(midx) @@ -140,7 +142,7 @@ def symmetry_index(d, r, dtype=None): #print(midx0-midx) #midx = midx0 - P = bm.concatenate([bm.tensor([1]), bm.cumprod(bm.arange(r+1)[1:], axis=0)], + P = bm.concatenate([bm.tensor([1],device=device), bm.cumprod(bm.arange(r+1, device=device)[1:], axis=0)], axis=0, dtype=dtype) num = P[r]/bm.prod(P[midx], axis=1) return symidx, num diff --git a/fealpy/mesh/triangle_mesh.py b/fealpy/mesh/triangle_mesh.py index 882661075..cd5360a7d 100644 --- a/fealpy/mesh/triangle_mesh.py +++ b/fealpy/mesh/triangle_mesh.py @@ -17,12 +17,13 @@ def __init__(self, node: TensorLike, cell: TensorLike) -> None: """ super().__init__(TD=2, itype=cell.dtype, ftype=node.dtype) kwargs = bm.context(cell) - self.device = device - + self.node = node self.cell = cell - self.localEdge = bm.tensor([(1, 2), (2, 0), (0, 1)], **kwargs, device=self.device) - self.localFace = bm.tensor([(1, 2), (2, 0), (0, 1)], **kwargs, device=self.device) + + + self.localEdge = bm.tensor([(1, 2), (2, 0), (0, 1)], **kwargs) + self.localFace = bm.tensor([(1, 2), (2, 0), (0, 1)], **kwargs) self.ccw = bm.tensor([0, 1, 2], **kwargs) self.localCell = bm.tensor([ @@ -1117,7 +1118,7 @@ def from_box(cls, box=[0, 1, 0, 1], nx=10, ny=10, *, threshold=None, itype = bm.int32 if ftype is None: ftype = bm.float64 - + NN = (nx + 1) * (ny + 1) x = bm.linspace(box[0], box[1], nx+1, dtype=ftype, device=device) y = bm.linspace(box[2], box[3], ny+1, dtype=ftype, device=device) diff --git a/fealpy/pde/biharmonic_triharmonic_2d.py b/fealpy/pde/biharmonic_triharmonic_2d.py index 846fed63e..432fda355 100644 --- a/fealpy/pde/biharmonic_triharmonic_2d.py +++ b/fealpy/pde/biharmonic_triharmonic_2d.py @@ -59,7 +59,8 @@ class DoubleLaplacePDE(): """ \Delta^2 u = f """ - def __init__(self, u): + def __init__(self, u, device=None): + self.device = device x = sp.symbols("x") y = sp.symbols("y") @@ -91,10 +92,13 @@ def init_mesh(self, n=1, meshtype='poly'): def source(self, p): x = p[..., 0] y = p[..., 1] + x_cpu = bm.to_numpy(x) + y_cpu = bm.to_numpy(y) #sin = bm.sin #pi = bm.pi #cos = bm.cos - return self.L2u(x, y) + print(bm.array(self.L2u(x_cpu, y_cpu),device=self.device)) + return bm.array(self.L2u(x_cpu, y_cpu),device=self.device) def solution(self, p): x = p[..., 0]