（三）、Fealpy 创建有限元空间

这一讲，我们来介绍 Fealpy 中的有限元空间。 Fealpy 中实现了几种常用的有限元空间， 比如拉格朗日有限元空间， 缩放多项式空间等。
首先我们要知道fealpy 都是基于重心坐标来计算的。 因为 fealpy 想实现任意次的空间， 故我们要约定一下多重指标的排列：

给定多重指标的排列之后，我们便可以构造 $p$ 次多项式

注意这里的 ${\lambda }_i$ 是重心坐标， 其与一般坐标关系为

我们从一个具体的例子来看：

右下角的图显示了我们多重指标的排列方式。 将它的值除以三(p=3)，便得到了对应的重心坐标。图二显示了重心坐标编号在三角形中的对应方式。另外，总的局部自由度为 $(p+1)(p+2)/2$, $p=3$时， 为10.
并且，由有限元的知识，如果右下图的编号有两个值为0，比如编号0，6，9，那么这个编号对用的点便在顶点 $\cdots \cdots$
我们展示一个三维的编号图：

计算Lagrange基函数

如我们前面所说， 计算基函数便是要求如下式子：

在 $m_{\mathbf{\alpha }}!$ 的计算中， 我们要用到 $1!,2!,\cdots ,p!$, 我们还要用到 $p{\lambda }_0-1,p{\lambda }_0-2,\cdots ,p{\lambda }_0-(p-1),\cdots$ 等， 所以我们先构造两个数组

向量 $P$ 和矩阵 $A$ 包含我们计算基函数 ${\varphi }_{\mathbf{\alpha }}$ 的所有元素， 我们接下来考虑如何实现 ${\prod }_{j_1=0}^t(p{\lambda }_1-j_1),t\le m_1-1.$ 我们将矩阵第一行乘到第二行上， 第一二行乘到第三行上， 前n行乘到第 n+1 行上， 便可以得到想要的累乘式子， 再结合 $1/m_{\mathbf{\alpha }}!$，我们可以得到矩阵B


这样我们便实现了基函数的计算， 下面我们展示一下计算过程中的部分关键代码：

与计算基函数相同，计算基函数的导数可以通过下式完成：

具体细节可参加 魏老师 bilibili 录课视频。

首先我们从 Functionspace 中引入 LagrangeFiniteElementSpace, 当然也要引入之前介绍的网格工厂（Meshfactory）:

import numpy as np
from fealpy.mesh import MeshFactory
from fealpy.functionspace import LagrangeFiniteElementSpace
import matplotlib.pyplot as plt

利用 meshfactory 中的 boxmesh2d 来生成一个网格 ：

mf = MeshFactory()
box = [0, 1, 0, 1]
mesh = mf.boxmesh2d(box, nx=4, ny=4, meshtype='tri')

接下来调用 LagrangeFiniteElementSpace 生成有限元空间：

space = LagrangeFiniteElementSpace(mesh, p=3)

这里的 $p$ 是多项式次数， $p=1$ 为线性元， $p=2$ 为二次元， $p>2$ 为高次元。
我们可以查看一下生成网格的自由度

ldof = space.number_of_local_dofs()
gdof = space.number_of_global_dofs()

print('ldof:', ldof)
print('gdof:', gdof)

我们将图中的关系全部画出来，有

fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)
mesh.find_node(axes, showindex=True)
mesh.find_cell(axes, showindex=True)

接下来， 我们给定重心坐标，求解对应的基函数与基函数梯度

# 重心坐标
bc = np.array([[1/3, 1/3, 1/3]], dtype=np.float64)
print(mesh.bc_to_point(bc))  # 转化到真实网格点 (NQ, NC, 2)
phi = space.basis(bc)
gphi = space.grad_basis(bc)
print('重心坐标:', bc.shape)
print('phi:', phi.shape)
print('gphi:', gphi.shape)

结果为