Merge pull request #1237 from gaotingyi/master

光滑元三维插值

example/experiment/fem/interpolation_fem.py

@@ -0,0 +1,108 @@
+import argparse
+import sympy as sp
+from fealpy.backend import backend_manager as bm
+from fealpy.functionspace.cm_conforming_fe_space3d import CmConformingFESpace3d
+from fealpy.pde.biharmonic_triharmonic_3d import DoubleLaplacePDE, get_flist
+from fealpy.mesh.tetrahedron_mesh import TetrahedronMesh
+from fealpy.decorator import barycentric
+from fealpy.utils import timer
+
+## 参数解析
+parser = argparse.ArgumentParser(description=
+        """
+        光滑元有限元方法插值
+        """)
+
+parser.add_argument('--degree',
+        default=11, type=int,
+        help='光滑有限元空间的次数, 默认为 11 次.')
+
+parser.add_argument('--m',
+        default=1, type=int,
+        help='C^m 光滑元')
+
+parser.add_argument('--n',
+        default=1, type=int,
+        help='初始网格剖分段数.')
+
+parser.add_argument('--maxit',
+        default=2, type=int,
+        help='默认网格加密求解的次数, 默认加密求解 4 次')
+
+parser.add_argument('--backend',
+        default='pytorch', type=str,
+        help='默认后端为numpy')
+
+args = parser.parse_args()
+tmr = timer()
+next(tmr)
+bm.set_backend(args.backend)
+decive = "cpu"
+p = args.degree
+m = args.m
+n = args.n
+maxit = args.maxit
+
+bm.set_default_device('cpu')
+x = sp.symbols('x')
+y = sp.symbols('y')
+z = sp.symbols('z')
+#u = sp.sin(2*x)*sp.sin(2*y)*sp.sin(z)
+u = (sp.sin(2*sp.pi*y)*sp.sin(2*sp.pi*x)*sp.sin(2*sp.pi*z))**4
+flist = get_flist(u)
+NDof = bm.zeros(maxit, dtype=bm.float64)
+
+errorType = ['$|| u - u_h||_{\\Omega,0}$',
+             '$||\\nabla u - \\nabla u_h||_{\\Omega,0}$',
+             '$||\\nabla^2 u - \\nabla^2 u_h||_{\\Omega,0}$',
+             '$||\\nabla^3 u - \\nabla^3 u_h||_{\\Omega,0}$']
+errorMatrix = bm.zeros((4, maxit), dtype=bm.float64)
+tmr.send('网格和pde生成时间')
+
+for i in range(maxit):
+    mesh = TetrahedronMesh.from_box([0,1,0,1,0,1], 2**i, 2**i, 2**i)
+    node = mesh.entity('node')
+    isCornerNode = bm.zeros(len(node),dtype=bm.bool)
+    for n in bm.array([[0,0,0],[1,0,0],[0,1,0],[1,1,0],[1,1,1],[0,0,1],[1,0,1],[0,1,1]], dtype=bm.float64):
+        isCornerNode = isCornerNode | (bm.linalg.norm(node-n[None, :], axis=1)<1e-10)
+    mesh.node = node/2
+    space = CmConformingFESpace3d(mesh, p=p, m=m, isCornerNode=isCornerNode)
+    tmr.send(f'第{i}次空间生成时间')
+    fI = space.interpolation(flist)
+    tmr.send(f'第{i}次插值时间')
+
+    @barycentric
+    def ug1val(p):
+        return fI.grad_m_value(p, 1)
+
+    @barycentric
+    def ug2val(p):
+        return fI.grad_m_value(p, 2)
+    @barycentric
+    def ug3val(p):
+        return fI.grad_m_value(p, 3)
+
+    errorMatrix[0, i] = mesh.error(flist[0], fI, q=p+3)
+    import ipdb
+    ipdb.set_trace()
+    errorMatrix[1, i] = mesh.error(flist[1], ug1val, q=p+3)
+    #errorMatrix[2, i] = mesh.error(flist[2], ug2val, q=p+3)
+    #errorMatrix[3, i] = mesh.error(flist[3], ug3val, q=p+3)
+    print(errorMatrix)
+
+
+next(tmr)
+
+
+next(tmr)
+print("最终误差",errorMatrix)
+print("order : ", bm.log2(errorMatrix[0,:-1]/errorMatrix[0,1:]))
+print("order : ", bm.log2(errorMatrix[1,:-1]/errorMatrix[1,1:]))
+print("order : ", bm.log2(errorMatrix[2,:-1]/errorMatrix[2,1:]))
+
+
+
+
+
+
+

fealpy/functionspace/cm_conforming_fe_space3d.py

@@ -20,12 +20,13 @@
 _S = slice(None)
 
 class CmConformingFESpace3d(FunctionSpace, Generic[_MT]):
-    def __init__(self, mesh:_MT, p: int, m: int, isCornerNode:bool):
+    def __init__(self, mesh:_MT, p: int, m: int, isCornerNode:bool): 
         assert(p>8*m)
         self.mesh = mesh
         self.p = p
         self.m = m
         self.bspace = BernsteinFESpace(mesh, p)
+        self.device = mesh.device
 
         self.ftype = mesh.ftype
         self.itype = mesh.itype

fealpy/pde/biharmonic_triharmonic_3d.py

@@ -1,4 +1,4 @@
-import numpy as np
+from fealpy.backend import backend_manager as bm
 import sympy as sp
 from fealpy.mesh import TetrahedronMesh 
 #from fealpy.mesh import TriangleMesh
@@ -34,9 +34,9 @@ def source(self, p):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        sin = np.sin
-        pi = np.pi
-        cos = np.cos
+        sin = bm.sin
+        pi = bm.pi
+        cos = bm.cos
         return -self.Lu(x, y, z)
 
     def solution(self, p):
@@ -49,10 +49,10 @@ def gradient(self, p):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        sin = np.sin
-        pi = np.pi
-        cos = np.cos
-        val = np.zeros(p.shape, dtype=np.float64)
+        sin = bm.sin
+        pi = bm.pi
+        cos = bm.cos
+        val = bm.zeros(p.shape, dtype=bm.float64)
         val[..., 0] = self.ux(x, y)
         val[..., 1] = self.uy(x, y) 
         val[..., 2] = self.uz(x, y)
@@ -116,10 +116,10 @@ def gradient(self, p):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        sin = np.sin
-        pi = np.pi
-        cos = np.cos
-        val = np.zeros(p.shape, dtype=np.float64)
+        sin = bm.sin
+        pi = bm.pi
+        cos = bm.cos
+        val = bm.zeros(p.shape, dtype=bm.float64)
         val[..., 0] = self.ux(x, y, z)
         val[..., 1] = self.uy(x, y, z) 
         val[..., 2] = self.uz(x, y, z)
@@ -130,10 +130,10 @@ def hessian(self, p):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        sin = np.sin
-        pi = np.pi
-        cos = np.cos
-        val = np.zeros(p.shape[:-1]+(6, ), dtype=np.float64)
+        sin = bm.sin
+        pi = bm.pi
+        cos = bm.cos
+        val = bm.zeros(p.shape[:-1]+(6, ), dtype=bm.float64)
         val[..., 0] = self.uxx(x, y, z) 
         val[..., 1] = self.uxy(x, y, z) 
         val[..., 2] = self.uxz(x, y, z)
@@ -216,9 +216,9 @@ def source(self, p):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        sin = np.sin
-        pi = np.pi
-        cos = np.cos
+        sin = bm.sin
+        pi = bm.pi
+        cos = bm.cos
         return -self.L3u(x, y, z) 
 
     def solution(self, p):
@@ -231,10 +231,10 @@ def gradient(self, p):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        sin = np.sin
-        pi = np.pi
-        cos = np.cos
-        val = np.zeros(p.shape, dtype=np.float64)
+        sin = bm.sin
+        pi = bm.pi
+        cos = bm.cos
+        val = bm.zeros(p.shape, dtype=bm.float64)
         val[..., 0] = self.ux(x, y, z)
         val[..., 1] = self.uy(x, y, z) 
         val[..., 2] = self.uz(x, y, z)
@@ -244,10 +244,10 @@ def hessian(self, p):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        sin = np.sin
-        pi = np.pi
-        cos = np.cos
-        val = np.zeros(p.shape[:-1]+(6, ), dtype=np.float64)
+        sin = bm.sin
+        pi = bm.pi
+        cos = bm.cos
+        val = bm.zeros(p.shape[:-1]+(6, ), dtype=bm.float64)
         val[..., 0] = self.uxx(x, y, z) 
         val[..., 1] = self.uxy(x, y, z) 
         val[..., 2] = self.uxz(x, y, z)
@@ -260,10 +260,10 @@ def grad_3(self, p):
         x = p[..., 0]
         y = p[..., 1]
         z = p[..., 2]
-        sin = np.sin
-        pi = np.pi
-        cos = np.cos
-        val = np.zeros(p.shape[:-1]+(10, ), dtype=np.float64)
+        sin = bm.sin
+        pi = bm.pi
+        cos = bm.cos
+        val = bm.zeros(p.shape[:-1]+(10, ), dtype=bm.float64)
 
         val[..., 0] = self.uxxx(x, y, z)              
         val[..., 1] = self.uxxy(x, y, z)
@@ -365,32 +365,35 @@ def get_flist(u_sp):
     uzzzz = sp.lambdify(('x', 'y', 'z'), uzzzz_sp, 'numpy')
 
     f    = lambda node : u(node[..., 0], node[..., 1], node[..., 2])
+
+    @cartesian
     def grad_f(node):
         x = node[..., 0]
         y = node[..., 1]
         z = node[..., 2]
-        val = np.zeros_like(node)
+        val = bm.zeros_like(node)
         val[..., 0] = ux(x, y, z)
         val[..., 1] = uy(x, y, z)
         val[..., 2] = uz(x, y, z)
-        return val 
+        return bm.array(val) 
+    @cartesian
     def grad_2_f(node):
         x = node[..., 0]
         y = node[..., 1]
         z = node[..., 2]
-        val = np.zeros(x.shape+(6, ), dtype=np.float64)
+        val = bm.zeros(x.shape+(6, ), dtype=bm.float64)
         val[..., 0] = uxx(x, y, z) 
         val[..., 1] = uxy(x, y, z) 
         val[..., 2] = uxz(x, y, z)
         val[..., 3] = uyy(x, y, z) 
         val[..., 4] = uyz(x, y, z) 
         val[..., 5] = uzz(x, y, z)
-        return val 
+        return bm.array(val) 
     def grad_3_f(node):
         x = node[..., 0]
         y = node[..., 1]
         z = node[..., 2]
-        val = np.zeros(x.shape+(10, ), dtype=np.float64)
+        val = bm.zeros(x.shape+(10, ), dtype=bm.float64)
         val[..., 0] = uxxx(x, y, z)              
         val[..., 1] = uxxy(x, y, z)
         val[..., 2] = uxxz(x, y, z) 
@@ -406,7 +409,7 @@ def grad_4_f(node):
         x = node[..., 0]
         y = node[..., 1]
         z = node[..., 2]
-        val = np.zeros(x.shape+(15, ), dtype=np.float64)
+        val = bm.zeros(x.shape+(15, ), dtype=bm.float64)
         val[..., 0] = uxxxx(x, y, z)
         val[..., 1] = uyxxx(x, y, z)
         val[..., 2] = uxxxz(x, y, z)