Merge branch 'master' of github.com:weihuayi/fealpy

fealpy/fem/LinearElasticityOperatorIntegrator.py

@@ -46,7 +46,10 @@ def assembly_cell_matrix(self, space, index=np.s_[:], cellmeasure=None, out=None
 
         A = [np.einsum('i, ijm, ijn, j->jmn', ws, grad[..., i], grad[..., j], cellmeasure, optimize=True) for i, j in idx]
 
-        D = mu*np.sum(A)
+        D = 0
+        for i in range(GD):
+            D += mu*A[imap[(i, i)]]
+
         if space[0].doforder == 'sdofs': # 标量自由度优先排序 
             for i in range(GD):
                 for j in range(i, GD):

fealpy/functionspace/FirstNedelecFiniteElementSpace2d.py

@@ -54,12 +54,12 @@ def is_boundary_dof(self, threshold=None):
         isBdDof[edge2dof[idx]] = True
         return isBdDof
 
-    def edge_to_dof(self):
+    def edge_to_dof(self, index=np.s_[:]):
         mesh = self.mesh
         NE = mesh.number_of_edges()
         edof = self.number_of_local_dofs('edge')
         edge2dof = np.arange(NE * edof).reshape(NE, edof)
-        return edge2dof
+        return edge2dof[index]
 
     def cell_to_dof(self):
         """
@@ -142,7 +142,7 @@ def basis(self, bc, index=np.s_[:]):
         -----
         """
 
-        # 每个单元上的全部自由度个数
+                # 每个单元上的全部自由度个数
         ldof = self.number_of_local_dofs(doftype='all')
         edof = self.number_of_local_dofs(doftype='edge')
 
@@ -359,8 +359,12 @@ def mass_matrix(self, c=None, q=None):
             if isinstance(c, (int, float)):
                 M = np.einsum('i, ijkd, ijmd, j->jkm', c * ws, phi, phi,
                               self.cellmeasure, optimize=True)
-            else:
-                M = np.einsum('i, ij, ijkd, ijmd, j->jkm', ws, c, phi, phi, cellmeasure, optimize=True)
+            elif isinstance(c, np.ndarray):
+                if len(c.shape)==2:
+                    M = np.einsum('i, ij, ijkd, ijmd, j->jkm', ws, c, phi, phi, cellmeasure, optimize=True)
+                elif len(c.shape)==4:
+                    print("aaa")
+                    M = np.einsum('i, ijdl, ijkd, ijml, j->jkm', ws, c, phi, phi, cellmeasure, optimize=True)
         cell2dof = self.cell_to_dof()
         gdof = self.number_of_global_dofs()
 
@@ -370,6 +374,46 @@ def mass_matrix(self, c=None, q=None):
         M = csr_matrix((M.flat, (I.flat, J.flat)), shape=(gdof, gdof))
         return M
 
+    def edge_basis(self, bcs, index=np.s_[:]):
+        """
+        @brief 计算每条边上的 v \cdot t(即 n \times v) 在重心坐标处的值.
+                对于线性元，这个值就是 1/e
+        """
+        em = self.mesh.entity_measure("edge", index=index)
+        shape = bcs.shape[:-1] + (len(em), )
+        val = np.broadcast_to(1/em, shape,) #(NQ, NE)
+        return val
+
+    def edge_mass_matrix(self, c = 1, index=np.s_[:]):
+        """
+        @brief (n \times u, n \times v)_{Gamma_{robin}}
+               注意 : n_i \times \phi_i = 1/e_i
+        @param c 系数, 现在只考虑了 c 是常数的情况
+        """
+        edge2dof = self.dof.edge_to_dof(index=index)
+        em = self.mesh.entity_measure("edge", index=index)
+        gdof = self.dof.number_of_global_dofs()
+
+        EM = c*csr_matrix((em.flat, (edge2dof.flat, edge2dof.flat)), shape=(gdof, gdof))
+        return EM
+
+    def robin_vector(self, f, isRobinEdge):
+        """
+        @brief 计算 (f, n\times v)_{\Gamma_{robin}} 其中 n \times v = 1/e
+        """
+        eint = self.integralalg.edgeintegrator
+        bcs, ws = eint.get_quadrature_points_and_weights()
+        n = self.mesh.edge_unit_normal()[isRobinEdge]
+
+        point = self.mesh.bc_to_point(bcs, index=isRobinEdge)
+        fval = f(point, n)
+
+        e2dof = self.dof.edge_to_dof(index=isRobinEdge)
+        gdof = self.dof.number_of_global_dofs()
+        F = np.zeros(gdof, dtype=np.float_)
+        F[e2dof] = np.einsum("qe, q->e", fval, ws)[:, None] 
+        return F
+
     def curl_matrix(self, c=None, q=None):
         """
 
@@ -399,6 +443,7 @@ def curl_matrix(self, c=None, q=None):
             if isinstance(c, (int, float)):
                 M = np.einsum('i, ijk, ijm, j->jkm', c * ws, phi, phi, cellmeasure, optimize=True)
             else:
+                print("aaa")
                 M = np.einsum('q, qc, qck..., qcm..., c->ckm', ws, c, phi, phi, cellmeasure, optimize=True)
 
         cell2dof = self.cell_to_dof()
@@ -416,6 +461,9 @@ def source_vector(self, f):
         b = self.integralalg.construct_vector_v_v(f, self.basis, cell2dof, gdof=gdof)
         return b
 
+    def face_mass_matrix(self):
+        pass
+
     def set_dirichlet_bc(self, gD, uh, threshold=None, q=None):
         """
         """

fealpy/functionspace/FirstNedelecFiniteElementSpace3d.py

@@ -322,12 +322,11 @@ def set_dirichlet_bc(self, gD, uh, threshold=None, q=None):
         if 1: #节点型自由度
             locEdge = np.array([[1, 2], [2, 0], [0, 1]], dtype=np.int_)
             point = 0.5*(np.sum(node[face[:, locEdge][index]], axis=-2))
-            gval = gD(point) #(NF, 3, 3)
-
-            vec = mesh.edge_tangent()[face2edge]#(NF, 3, 3)
+            vec = mesh.edge_tangent()[face2edge]
+            gval = gD(point, vec) #(NF, 3)
 
             face2dof = self.dof.face_to_dof()[index]
-            uh[face2dof] = np.sum(gval*vec, axis=-1) 
+            uh[face2dof] = gval 
         else: #积分型自由度
             bcs, ws = self.integralalg.edgeintegrator.get_quadrature_points_and_weights()
             ps = mesh.bc_to_point(bcs)[:, face2edge]

fealpy/geometry/implicit_curve.py

@@ -196,10 +196,6 @@ def value(self, p):
         return self(p)
 
     def gradient(self, p):
-        k = self.k
-        b = self.b
-        l = self.l
-
         val = self(p)
 
         x = p[..., 0] - self.center[0]

fealpy/mesh/HalfEdgeMesh2d.py

@@ -465,6 +465,9 @@ def entity(self, etype=2):
         else:
             raise ValueError("`entitytype` is wrong!")
 
+    def geo_dimension(self):
+        return self.node.shape[1]
+
     def entity_barycenter(self, etype='cell', index=np.s_[:]):
         node = self.entity('node')
         GD = self.geo_dimension()

fealpy/mesh/StructureMesh3dDataStructure.py

@@ -124,54 +124,50 @@ def face2cell(self):
 
         # x direction
         NF0 = 0
-        NF1 = ny * nz
-        face2cell[NF0:NF1, 0] = idx[0].flatten()
-        face2cell[NF0:NF1, 1] = idx[0].flatten()
-        face2cell[NF0:NF1, 2:4] = 0
+        NF1 = (nx+1) * ny * nz
+        face2cell[NF0:NF1-ny*nz, 0] = idx.flatten()
+        face2cell[NF0+ny*nz:NF1, 1] = idx.flatten()
+        face2cell[NF0:NF1-ny*nz, 2] = 0
+        face2cell[NF0:NF1-ny*nz, 3] = 1
 
-        NF0 = NF1
-        NF1 += nx * ny * nz
-        face2cell[NF0:NF1, 0] = idx.flatten()
-        face2cell[NF0:NF1, 2] = 1
-        face2cell[NF0:NF1 - ny * nz, 1] = idx[1:].flatten()
-        face2cell[NF0:NF1 - ny * nz, 3] = 0
-        face2cell[NF1 - ny * nz:NF1, 1] = idx[-1].flatten()
-        face2cell[NF1 - ny * nz:NF1, 3] = 1
+        face2cell[NF1-ny*nz:NF1, 0] = idx[-1].flatten()
+        face2cell[NF0:NF0+ny*nz, 1] = idx[0].flatten()
+        face2cell[NF1-ny*nz:NF1, 2] = 1
+        face2cell[NF0:NF0+ny*nz, 3] = 0
 
         # y direction
-        c = np.transpose(idx, (1, 0, 2))
+        idy = np.swapaxes(idx, 1, 0)
         NF0 = NF1
-        NF1 += nx * nz
-        face2cell[NF0:NF1, 0] = c[0].flatten()
-        face2cell[NF0:NF1, 1] = c[0].flatten()
-        face2cell[NF0:NF1, 2:4] = 2
+        NF1 += nx * (ny+1) * nz
 
-        NF0 = NF1
-        NF1 += nx * ny * nz
-        face2cell[NF0:NF1, 0] = c.flatten()
-        face2cell[NF0:NF1, 2] = 3
-        face2cell[NF0:NF1 - nx * nz, 1] = c[1:].flatten()
-        face2cell[NF0:NF1 - nx * nz, 3] = 2
-        face2cell[NF1 - nx * nz:NF1, 1] = c[-1].flatten()
-        face2cell[NF1 - nx * nz:NF1, 3] = 3
+        fidy = np.arange(NF0, NF1).reshape(nx, ny+1, nz).swapaxes(0, 1)
+
+        face2cell[fidy[:-1], 0] = idy
+        face2cell[fidy[1:], 1] = idy
+        face2cell[fidy[:-1], 2] = 0
+        face2cell[fidy[1:], 3] = 1
+
+        face2cell[fidy[-1], 0] = idy[-1]
+        face2cell[fidy[0], 1] = idy[0]
+        face2cell[fidy[-1], 2] = 1
+        face2cell[fidy[0], 3] = 0
 
         # z direction
-        c = np.transpose(idx, (2, 0, 1))
+        idz = np.transpose(idx, (2, 0, 1))
         NF0 = NF1
-        NF1 += nx * ny
-        face2cell[NF0:NF1, 0] = c[0].flatten()
-        face2cell[NF0:NF1, 1] = c[0].flatten()
-        face2cell[NF0:NF1, 2:4] = 4
+        NF1 += nx * ny * (nz + 1)
 
-        NF0 = NF1
-        NF1 += nx * ny * nz
-        face2cell[NF0:NF1, 0] = c.flatten()
-        face2cell[NF0:NF1, 2] = 5
-        face2cell[NF0:NF1 - nx * ny, 1] = c[1:].flatten()
-        face2cell[NF0:NF1 - nx * ny, 3] = 4
-        face2cell[NF1 - nx * ny:NF1, 1] = c[-1].flatten()
-        face2cell[NF1 - nx * ny:NF1, 3] = 5
+        fidz = np.arange(NF0, NF1).reshape(nx, ny, nz+1).transpose(2, 0, 1)
+
+        face2cell[fidz[:-1], 0] = idz
+        face2cell[fidz[1:], 1] = idz
+        face2cell[fidz[:-1], 2] = 0
+        face2cell[fidz[1:], 3] = 1
 
+        face2cell[fidz[-1], 0] = idz[-1]
+        face2cell[fidz[0], 1] = idz[0]
+        face2cell[fidz[-1], 2] = 1
+        face2cell[fidz[0], 3] = 0
         return face2cell
 
     @property

fealpy/mesh/UniformMesh3d.py

@@ -659,6 +659,16 @@ def entity(self, etype):
 
         @throws ValueError if the given etype is invalid.
         """
+        if etype in {'cell', 3}:
+            return self.ds.cell
+        elif etype in {'face', 2}:
+            return self.ds.face
+        elif etype in {'edge', 1}:
+            return self.ds.edge
+        elif etype in {'node', 0}:
+            return self.node.reshape(-1, 3)
+        else:
+            raise ValueError("`etype` is wrong!")
         pass
 
     ## @ingroup FEMInterface

fealpy/pde/MaxwellPDE2d.py

@@ -5,7 +5,13 @@
 from sympy.vector import CoordSys3D, Del, curl
 
 class MaxwellPDE2d():
-    def __init__(self, f):
+    def __init__(self, f, beta = 1, k = 1):
+        """
+        @brief 求解方程 
+                curl curl E - beta E = J     Omega
+                          n \times E = g     Gamma0
+             - curl E - k n \times E = f     Gamma1
+        """
         C = CoordSys3D('C')
         x = sym.symbols("x")
         y = sym.symbols("y")
@@ -29,6 +35,9 @@ def __init__(self, f):
         self.curlcurlFx = sym.lambdify(('x', 'y'), ccfx, "numpy")
         self.curlcurlFy = sym.lambdify(('x', 'y'), ccfy, "numpy")
 
+        self.beta = beta
+        self.k = k
+
     @cartesian
     def solution(self, p):
         x = p[..., 0, None]
@@ -75,7 +84,7 @@ def source(self, p):
         if type(ccFy) is not np.ndarray:
             ccFy = np.ones(x.shape, dtype=np.float_)*ccFy
         ccf = np.c_[ccFx, ccFy] 
-        return ccf - self.solution(p)
+        return ccf - self.beta * self.solution(p)
 
     @cartesian
     def dirichlet(self, p, t):
@@ -88,31 +97,35 @@ def neumann(self, p, n):
         @param p: (..., N, ldof, 3)
         @param n: (N, 3)
         """
-        x = p[..., 0, None]
-        y = p[..., 1, None]
-        cFx = self.curlFx(x, y)
-        cFy = self.curlFy(x, y)
-        if type(cFx) is not np.ndarray:
-            cFx = np.ones(x.shape, dtype=np.float_)*cFx
-        if type(cFy) is not np.ndarray:
-            cFy = np.ones(x.shape, dtype=np.float_)*cFy
-        cf = np.c_[cFx, cFy] #(..., NC, ldof, 3)
-        return np.cross(n[None, :], cf)
+        pass
+
+    def robin(self, p, n):
+        """!
+        @param p: (..., N, ldof, 3)
+        @param n: (N, 3)
+        """
+        cf = self.curl_solution(p)
+        fval = self.solution(p)
+        return -cf - np.cross(n[None, :], fval)
 
     def boundary_type(self, mesh):
         bdface = mesh.ds.boundary_face_index()
         f2n = mesh.face_unit_normal()[bdface]
-        neumannbd = np.abs(f2n[:, 2])<0.9
-        bd = {"neumann": bdface[neumannbd], "dirichlet": bdface[~neumannbd]}
+        isLeftBd   = np.abs(f2n[:, 0]+1)<1e-14
+        isRightBd  = np.abs(f2n[:, 0]-1)<1e-14
+        isBottomBd = np.abs(f2n[:, 1]+1)<1e-14
+        isUpBd     = np.abs(f2n[:, 1]-1)<1e-14
+        bd = {"left": bdface[isLeftBd], "right": bdface[isRightBd], 
+              "up": bdface[isUpBd], "bottom": bdface[isBottomBd], }
         return bd
 
 class SinData(MaxwellPDE2d):
-    def __init__(self):
+    def __init__(self, beta = 1, k = 1):
         C = CoordSys3D('C')
         f = sym.sin(sym.pi*C.y)*C.i + sym.sin(sym.pi*C.x)*C.j
         #f = sym.sin(C.y)*C.i + sym.sin(C.x)*C.j
         #f = C.x*C.y*(1-C.x)*(1-C.y)*C.i + sym.sin(sym.pi*C.x)*sym.sin(sym.pi*C.y)*C.j
-        super(SinData, self).__init__(f)
+        super(SinData, self).__init__(f, beta, k)
 
     def init_mesh(self, nx=1, ny=1, meshtype='tri'):
         box = [0, 1, 0, 1]

test/fem/test_bilinear_form.py

@@ -59,7 +59,7 @@ def test_tetrahedron_mesh():
 
     p=1
     mesh = TetrahedronMesh.from_one_tetrahedron()
-    mesh.uniform_refine()
+    mesh.uniform_refine(n=1)
     space = Space(mesh, p=p)
 
     bform = BilinearForm(space)
@@ -127,25 +127,29 @@ def test_linear_elasticity_model():
     ny = 10 
     p = 1 
 
-    mesh = TriangleMesh.from_unit_square(nx=nx*2, ny=ny*2)
+#    mesh = TriangleMesh.from_unit_square(nx=nx*2, ny=ny*2)
+    mesh = TriangleMesh.from_one_triangle()
     space = Space(mesh, p=p, doforder='sdofs')
     ospace = OldSpace(mesh, p=p)
 
     GD = mesh.geo_dimension()
     bform = BilinearForm(GD*(space,))
     bform.add_domain_integrator(LinearElasticityOperatorIntegrator(lam=pde.lam,
-        mu=pde.mu, q=p+2))
+        mu=pde.mu, q=p))
 
     bform.assembly()
 
     A = bform.get_matrix()
-    B = ospace.linear_elasticity_matrix(pde.lam, pde.mu, q=p+2)
+    B = ospace.linear_elasticity_matrix(pde.lam, pde.mu, q=p)
+    print(A.toarray())
+    print(B.toarray())
+
     np.testing.assert_array_almost_equal(A.toarray(), B.toarray())
 
 
 if __name__ == '__main__':
-    test_linear_elasticity_model()
-    #test_tetrahedron_mesh()
+    #test_linear_elasticity_model()
+    test_tetrahedron_mesh()
     #test_triangle_mesh()
     #test_truss_structure()