Merge pull request #1220 from july-liuzhi/master

update

example/experiment/fem/maxwell_eq_n1fem_dirichlet.py

@@ -64,7 +64,9 @@
 args = parser.parse_args()
 p = args.degree
 maxit = args.maxit
+maxit = 3
 dim = args.dim
+dim = 3
 backend = args.backend
 bm.set_backend(backend)
 if dim == 2:

example/experiment/fem/maxwell_eq_n2fem_dirichlet.py

@@ -65,9 +65,11 @@
 args = parser.parse_args()
 p = args.degree
 maxit = args.maxit
+maxit = 3
 dim = args.dim
+dim = 3
 backend = args.backend
-bm.set_backend("pytorch")
+bm.set_backend(backend)
 if dim == 2:
     pde = PDE2d()
 else:
@@ -79,8 +81,8 @@
              '$||\\nabla \\times u - \\nabla_h \\times u_h||_0$ with k=2',
              '$||\\nabla \\times u - \\nabla_h \\times u_h||_0$ with k=3',
              '$||\\nabla \\times u - \\nabla_h \\times u_h||_0$ with k=4']
-errorMatrix = bm.zeros((len(errorType), maxit), dtype=bm.float64)
-NDof = bm.zeros(maxit, dtype=bm.float64)
+errorMatrix = np.zeros((len(errorType), maxit), dtype=np.float64)
+NDof = np.zeros(maxit, dtype=np.float64)
 
 tmr = timer()
 next(tmr)

fealpy/functionspace/first_nedelec_fe_space_3d.py

@@ -24,6 +24,7 @@ def __init__(self, mesh, p):
         self.multiindex = mesh.multi_index_matrix(p,3)
         self.ftype = mesh.ftype
         self.itype = mesh.itype
+        self.device = mesh.device
 
     def number_of_local_dofs(self, doftype='all'):
         p = self.p
@@ -74,7 +75,7 @@ def face_to_dof(self):
         edge = self.mesh.entity('edge')
         face = self.mesh.entity('face')
 
-        f2d = bm.zeros((NF, fdof), dtype=self.itype)
+        f2d = bm.zeros((NF, fdof), device=self.device,dtype=self.itype)
         #f2d[:, :eldof*3] = e2dof[f2e].reshape(NF, eldof*3)
         f2d = bm.set_at(f2d,(slice(None),slice(eldof*3)),e2dof[f2e].reshape(NF, eldof*3))
         s = [1, 0, 0]
@@ -224,7 +225,7 @@ def basis(self, bcs, index=_S):
 
         c2esign = mesh.cell_to_edge_sign() #(NC, 3, 2)
 
-        l = bm.zeros((4, )+bcs[None,: ,0,None, None].shape, dtype=self.ftype)
+        l = bm.zeros((4, )+bcs[None,: ,0,None, None].shape,device=self.device, dtype=self.ftype)
 
         l = bm.set_at(l,(0),bcs[None, :,0,None, None])
         l = bm.set_at(l,(1),bcs[None, :,1,None, None])
@@ -300,7 +301,7 @@ def face_internal_basis(self, bcs, index=_S):
 
         fdof = self.dof.number_of_local_dofs("face")
         glambda = mesh.grad_face_lambda()
-        val = bm.zeros((NF,) + bcs.shape[:-1]+(fdof, 3), dtype=self.ftype)
+        val = bm.zeros((NF,) + bcs.shape[:-1]+(fdof, 3),device=self.device, dtype=self.ftype)
 
         l = bm.zeros((3, )+bcs[None, :,0, None, None].shape, dtype=self.ftype)
         l = bm.set_at(l,(0),bcs[None, :,0,  None, None])
@@ -338,7 +339,7 @@ def face_basis(self, bcs, index=_S):
         # edge basis
         phi = self.bspace.basis(bcs, p=p)
         multiIndex = self.mesh.multi_index_matrix(p, 2)
-        val = bm.zeros((NF,) + bcs.shape[:-1]+(ldof, 3), dtype=self.ftype)
+        val = bm.zeros((NF,) + bcs.shape[:-1]+(ldof, 3),device=self.device, dtype=self.ftype)
         for i in range(3):
             phie = phi[:, :, multiIndex[:, i]==0]
             c2esi = c2esign[:, i]
@@ -384,7 +385,7 @@ def curl_basis(self, bcs):
         phi = self.bspace.basis(bcs, p=p)
         gphi = self.bspace.grad_basis(bcs, p=p)
         multiIndex = self.mesh.multi_index_matrix(p, 3)
-        val = bm.zeros((NC,)+bcs.shape[:-1]+(ldof, 3), dtype=self.ftype)
+        val = bm.zeros((NC,)+bcs.shape[:-1]+(ldof, 3),device=self.device, dtype=self.ftype)
 
         # edge basis
         locEdgeDual = bm.tensor([[2, 3], [1, 3], [1, 2], [0, 3], [0, 2], [0, 1]])

fealpy/functionspace/second_nedelec_fe_space_2d.py

@@ -96,7 +96,7 @@ def cell_to_dof(self):
         c2e = self.mesh.cell_to_edge()
         c2esign = self.mesh.cell_to_face_sign()
 
-        c2d = bm.zeros((NC, ldof), dtype=self.itype)
+        c2d = bm.zeros((NC, ldof), dtype=bm.int64)
         # c2d[:, e2ldof] : (NC, 3, p+1), e2dof[c2e] : (NC, 3, p+1)
         tmp = e2dof[c2e]
         # tmp[~c2esign] = tmp[~c2esign, ::-1]       

fealpy/functionspace/second_nedelec_fe_space_3d.py

@@ -26,7 +26,7 @@ def __init__(self, mesh, p):
         self.multiindex3 = mesh.multi_index_matrix(p,3)
         self.ftype = mesh.ftype
         self.itype = mesh.itype
-
+        self.device=mesh.device
 
     def edge_to_local_face_dof(self):
         multiindex = self.multiindex2
@@ -53,7 +53,7 @@ def face_to_local_dof(self):
         ldof = self.number_of_local_dofs()
 
         fdof = self.number_of_local_dofs('faceall')
-        f2ld = bm.zeros((4, fdof), dtype=self.itype)
+        f2ld = bm.zeros((4, fdof),device=self.device, dtype=self.itype)
         eldof = self.edge_to_local_face_dof()
         nldof = bm.tensor([[eldof[(i+1)%3, 0], eldof[(i+2)%3, -1]] for i in range(3)])
 
@@ -91,10 +91,10 @@ def cell_to_dof(self):
         cdof = self.number_of_local_dofs('cell')
         fdof = self.number_of_local_dofs('faceall')
 
-        isndof = bm.zeros(ldof, dtype=bm.bool)
+        isndof = bm.zeros(ldof, device=self.device, dtype=bm.bool)
         isndof[f2ld] = True
 
-        c2d = bm.zeros((NC, ldof), dtype=bm.int64)#
+        c2d = bm.zeros((NC, ldof), device=self.device, dtype=bm.int64)#
         idx = bm.zeros((NC, 3), dtype=self.itype)
         fe = bm.tensor([[0, 1], [0, 2], [1, 2]], dtype=self.itype) #局部边
         for i in range(4):
@@ -235,7 +235,7 @@ def basis(self, bc,index=_S):
         c2v = self.basis_vector() #(NC, ldof, GD)
 
         shape = bc.shape[:-1]
-        val = bm.zeros((NC,)+ shape + (ldof, GD), dtype=self.ftype)
+        val = bm.zeros((NC,)+ shape + (ldof, GD),device=self.device, dtype=self.ftype)
 
         bval = self.lspace.basis(bc) #(NC, NQ, ldof//3)
 
@@ -260,7 +260,7 @@ def face_basis(self, bc, index=_S):
         f2v = self.face_basis_vector(index=index)#(NF, ldof, GD)
         NF = len(f2v)
         shape = bc.shape[:-1]
-        val = bm.zeros((NF,) + shape+(ldof, GD), dtype=self.ftype)
+        val = bm.zeros((NF,) + shape+(ldof, GD), device=self.device, dtype=self.ftype)
 
         bval = self.lspace.basis(bc) #(NF, NQ, ldof//3)
 
@@ -293,7 +293,7 @@ def face_basis_vector(self, index=_S):
 
         NF = len(f2e)
 
-        f2v = bm.zeros((NF, ldof, GD), dtype=self.ftype)
+        f2v = bm.zeros((NF, ldof, GD),device=self.device, dtype=self.ftype)
         # f2v[:, :ldof//2] = e2t[:, 0, None] #(NF, ldof//2, 3)
         # f2v[:, ldof//2:] = e2n[:, 0, None]
         f2v = bm.set_at(f2v,(slice(None),slice(None,ldof//2)),e2t[:, 0, None])
@@ -340,7 +340,7 @@ def face_dof_vector(self, index=_S):
 
         NF = len(f2e)
 
-        f2v = bm.zeros((NF, ldof, GD), dtype=self.ftype)
+        f2v = bm.zeros((NF, ldof, GD), device=self.device,dtype=self.ftype)
         # f2v[:, :ldof//2] = e2t[:, 0, None] #(NF, ldof//2, 3)
         # f2v[:, ldof//2:] = e2n[:, 0, None]
         f2v = bm.set_at(f2v,(slice(None),slice(None,ldof//2)),e2t[:, 0, None])
@@ -366,7 +366,7 @@ def curl_basis(self, bc):
 
         c2v = self.basis_vector()#(NC, ldof, GD)
         sgval = self.lspace.grad_basis(bc) #(NQ, NC, lldof, GD)
-        val = bm.zeros((NC,)+(bc.shape[0], )+(ldof, GD), dtype=self.ftype)
+        val = bm.zeros((NC,)+(bc.shape[0], )+(ldof, GD),device=self.device, dtype=self.ftype)
 
         # val[..., :ldof//3, :] = bm.cross(sgval, c2v[:,None, :ldof//3, :])
         # val[..., ldof//3:2*(ldof//3), :] = bm.cross(sgval, c2v[:,None, ldof//3:2*(ldof//3), :])
@@ -544,7 +544,7 @@ def set_dirichlet_bc(self, gD, uh, threshold=None, q=None):
         uh[face2dof[:, ldof//2:]] = bm.sum(gval*f2v[:, ldof//2:], axis=-1)
         # uh[face2dof[:, :ldof//2]] =0
         # uh[face2dof[:, ldof//2:]] =0
-        isDDof = bm.zeros(gdof, dtype=bm.bool)
+        isDDof = bm.zeros(gdof, device=self.device,dtype=bm.bool)
         isDDof[face2dof] = True
         return uh,isDDof
 
@@ -564,7 +564,7 @@ def set_neumann_bc(self,gD):
         points = mesh.bc_to_point(bcs)[isbdFace]
         n = mesh.face_unit_normal()[isbdFace]
         hval = gD(points,n)
-        vec = bm.zeros(gdof, dtype=self.ftype)
+        vec = bm.zeros(gdof, device=self.device,dtype=self.ftype)
         k = bm.einsum('fqg, fqlg,q,f->fl', hval, bphi,ws,fm) # (NF, ldof)
         bm.add.at(vec,face2dof,k)
         #bm.scatter_add(vec,face2dof,k)