Merge pull request #1467 from july-liuzhi/develop

update Nedelec

fealpy/fem/face_mass_integrator.py

@@ -42,7 +42,7 @@ def fetch(self, space: _FS):
         q = space.p+3 if self.q is None else self.q
         qf = mesh.quadrature_formula(q, 'face')
         bcs, ws = qf.get_quadrature_points_and_weights()
-        phi = space.face_basis(bcs)
+        phi = space.face_basis(bcs,index=index)
         return bcs, ws, phi, facemeasure, index
 
     def assembly(self, space: _FS):

fealpy/functionspace/first_nedelec_fe_space_3d.py

@@ -60,7 +60,7 @@ def face_to_internal_dof(self, index=_S):
         fdof = self.number_of_local_dofs(doftype='face')
         return bm.arange(NE*edof, NE*edof+NF*fdof,device=self.device).reshape(NF, fdof)[index]
 
-    def face_to_dof(self):
+    def face_to_dof(self,index=_S):
         p = self.p
         fldof = self.number_of_local_dofs('face')
         eldof = self.number_of_local_dofs('edge')
@@ -87,7 +87,7 @@ def face_to_dof(self):
             #f2d = bm.set_at(f2d,(flag,slice(eldof*i,eldof*(i+1))),f2d[flag, eldof*i:eldof*(i+1)][:, ::-1])
         #f2d[:, eldof*3:] = self.face_to_internal_dof() 
         f2d = bm.set_at(f2d,(slice(None),slice(eldof*3,eldof*3+fldof)),self.face_to_internal_dof())
-        return f2d
+        return f2d[index]
 
     def cell_to_dof(self):
         p = self.p
@@ -466,6 +466,9 @@ def is_boundary_dof(self, threshold=None,method=None):
 
     def cell_to_dof(self):
         return self.dof.cell2dof
+
+    def face_to_dof(self,index=_S):
+        return self.dof.face_to_dof()[index]
 
     def number_of_global_dofs(self):
         return self.dof.number_of_global_dofs()
@@ -635,9 +638,11 @@ def set_dirichlet_bc(self, gd, uh, threshold=None, q=None,method=None):
             points = mesh.bc_to_point(bcs)[index1]
             n = mesh.face_unit_normal()[index1]
             n = n[:,None,:]
-            h2 = gd(points)
-            g = bm.cross(n, h2) 
-            g = bm.cross(g,n)
+            h2 = gd(points,n)
+
+            #g = bm.cross(n, h2) 
+            #g = bm.cross(g,n)
+            g = bm.cross(n,h2)
 
             g1 = bm.einsum("cqld, cqd,q,c->cl", fbasis, g, ws, fm)
 
@@ -663,7 +668,15 @@ def set_dirichlet_bc(self, gd, uh, threshold=None, q=None,method=None):
         Minv = Minv*em[:,None,None]
 
         points1 = mesh.bc_to_point(bcs)[index2]
-        h1 = gd(points1)
+
+        n1 = mesh.face_unit_normal()
+        n2 = bm.zeros((NE,3))
+        f2e = mesh.face_to_edge()
+        n2[f2e] = n1[:,None,:]
+        n2 = n2[:,None,:][index2]
+
+        h1 = gd(points1,n2)
+        h1 = bm.cross(n2,h1)
         b = bm.einsum('eqd, ed->eq', h1, t) 
 
         g2 = bm.einsum('eql, eq,q->el', bphi, b,ws)

fealpy/functionspace/second_nedelec_fe_space_3d.py

@@ -537,7 +537,10 @@ def set_dirichlet_bc(self, gd, uh, threshold=None, q=None,method=None):
         # gval = gD(point, nor[:, None]) #(NF, ldof//2, 3)
         # gval = bm.cross(gval, nor[:, None])
         # print(bm.max(bm.abs(gval)))
-        gval = gd(point)
+        n = mesh.face_unit_normal()[index]
+        n = n[:,None,:]
+        gval = gd(point,n)
+        gval = bm.cross(n,gval)
 
 
         uh[face2dof[:, :ldof//2]] = bm.sum(gval*f2v[:, :ldof//2], axis=-1)
@@ -573,6 +576,9 @@ def set_neumann_bc(self,gD):
 
     def cell_to_dof(self):
         return self.dof.cell2dof
+
+    def face_to_dof(self):
+        return self.dof.face_to_dof()
 
     def number_of_global_dofs(self):
         return self.dof.number_of_global_dofs()

fealpy/pde/maxwell_3d.py

@@ -231,8 +231,9 @@ def source(self, p):
         return bm.concatenate([f1, f2, f3], axis=-1)
 
     @cartesian
-    def dirichlet(self, p):
-        val = self.solution(p)
+    def dirichlet(self, p,n):
+        val1 = self.solution(p)
+        val = bm.cross(val1,n)
         return val
 
     @cartesian