(七)、Fealpy 中的 PDE接口与实例

这一节， 我们来介绍一下 Fealpy 中的PDE接口。 打开 fealpy/fealpy/pde, 你可以看到大量的pde接口：

fealpy/example 中有很多求解的例子， 基本涵盖了 pde solver 等的东西，如下图：

因为作者的文档写的很清楚了， 故我这里只举一个例子来简单介绍一下：混合边界条件的poisson 方程求解
混合边界条件的方程如下：

# 构造一个二维有真解的例子，用 p 次有限元来求解， 画出误差阶， 输出误差表格
# -\Delta u + 3 u = f, in \Omega;
# u = gD  on \Gamma_D
# \partial u / \partial n = gN on  \Gamma_N
# \partial u / \partial n + u = gR on \Gamma_R
# 求解区域为 [0, 1]^2, 其中 $x = 0$ 是 Robin 边界条件;  x=1 是Neumann 边界条件, 其余为Dirichlet 边界条件,
# u(x,y) = x^2 + y^2 + cos(pi*x)cos(pi*y)

import numpy as np
from scipy.sparse.linalg import spsolve
from fealpy.mesh import MeshFactory
from fealpy.functionspace import LagrangeFiniteElementSpace
from fealpy.boundarycondition import DirichletBC, NeumannBC, RobinBC
from fealpy.decorator import cartesian, barycentric
from fealpy.tools.show import showmultirate, show_error_table
import matplotlib.pyplot as plt
class Model:

    def __init__(self):
        pass

    def domain(self):
        return np.array([0, 1, 0, 1])

    @cartesian
    def solution(self, p):
        """
		The exact solution
        Parameters
        ---------
        p :

        Examples
        -------
        p = np.array([0, 1], dtype=np.float64)
        p = np.array([[0, 1], [0.5, 0.5]], dtype=np.float64)
        """
        x = p[..., 0]
        y = p[..., 1]
        pi = np.pi
        val = np.cos(pi * x) * np.cos(pi * y)
        return val  # val.shape == x.shape

    @cartesian
    def source(self, p):
        """
		The right hand side of convection-diffusion-reaction equation
        INPUT:
            p: array object,
        """
        x = p[..., 0]
        y = p[..., 1]
        pi = np.pi
        val = (2 * pi ** 2 + 3) * np.cos(pi * x) * np.cos(pi * y)
        return val

    @cartesian
    def gradient(self, p):
        """
		The gradient of the exact solution
        """
        x = p[..., 0]
        y = p[..., 1]
        pi = np.pi
        val = np.zeros(p.shape, dtype=np.float64)
        val[..., 0] = -pi * np.sin(pi * x) * np.cos(pi * y)
        val[..., 1] = -pi * np.cos(pi * x) * np.sin(pi * y)
        return val  # val.shape == p.shape

    @cartesian
    def reaction_coefficient(self, p):
        x = p[..., 0]
        y = p[..., 1]
        return 3

    @cartesian
    def dirichlet(self, p):
        return self.solution(p)

    @cartesian
    def neumann(self, p, n):
        """
        Neuman  boundary condition

        Parameters
        ----------

        p: (NQ, NE, 2)
        n: (NE, 2)

        grad*n : (NQ, NE, 2)
        """
        grad = self.gradient(p)  # (NQ, NE, 2)
        val = np.sum(grad*n, axis=-1)  # (NQ, NE)
        return val

    @cartesian
    def robin(self, p, n):
        grad = self.gradient(p)  # (NQ, NE, 2)
        val = np.sum(grad*n, axis=-1)
        shape = len(val.shape)*(1, )
        kappa = np.array([1.0], dtype=np.float64).reshape(shape)
        val += self.solution(p)
        return val, kappa

    @cartesian
    def is_dirichlet_boundary(self, p):
        y = p[..., 1]
        return ( y == 1.0) | ( y == 0.0)

    @cartesian
    def is_neumann_boundary(self, p):
        x = p[..., 0]
        return x == 1.0

    @cartesian
    def is_robin_boundary(self, p):
        x = p[..., 0]
        return x == 0.0


# p = int(sys.argv[1])  # Lagrange 有限元多项式次数
# n = int(sys.argv[2])  # 初始网格剖分段数
# maxit = int(sys.argv[3])  # 网格加密最大次数

p = 2
n = 10
maxit = 4
pde = Model()
domain = pde.domain()

mf = MeshFactory()
mesh = mf.boxmesh2d(domain, nx=n, ny=n, meshtype='tri')

NDof = np.zeros(maxit, dtype=mesh.itype)
errorMatrix = np.zeros((2, maxit), dtype=mesh.ftype)
errorType = ['$|| u  - u_h ||_0$', '$|| \\nabla u - \\nabla u_h||_0$']
for i in range(maxit):
    print('Step:', i)
    space = LagrangeFiniteElementSpace(mesh, p=p)
    NDof[i] = space.number_of_global_dofs()
    uh = space.function()  # 返回一个有限元函数，初始自由度值全为 0
    A = space.stiff_matrix()
    M = space.mass_matrix(c=pde.reaction_coefficient)
    F = space.source_vector(pde.source)
    A += M
    bc1 = NeumannBC(space, pde.neumann, threshold=pde.is_neumann_boundary)
    F = bc1.apply(F=F)
    bc2 = RobinBC(space, pde.robin, threshold=pde.is_robin_boundary)
    A, F = bc2.apply(A=A, F=F)
    bc = DirichletBC(space, pde.dirichlet, threshold=pde.is_dirichlet_boundary)
    A, F = bc.apply(A, F, uh)

    uh[:] = spsolve(A, F)

    errorMatrix[0, i] = space.integralalg.error(pde.solution, uh.value, power=2)
    errorMatrix[1, i] = space.integralalg.error(pde.gradient, uh.grad_value,
                                                power=2)

    if i < maxit - 1:
        mesh.uniform_refine()

# 函数解图像
uh.add_plot(plt, cmap='rainbow')

# 收敛阶图像
showmultirate(plt, 0, NDof, errorMatrix, errorType, propsize=10)

# 输出误差的 latex 表格
show_error_table(NDof, errorType, errorMatrix)
plt.show()

结果为