(一)、Fealpy 创建各种各样的网格

Fealpy 库中实现了一个类 mesh 用于网格生成, 我们先介绍规则的网格（默认大家已经完成 fealpy 的安装， 见前言）
首先引入一些库：

import numpy as np
import matplotlib.pyplot as plt
from fealpy.mesh import QuadrangleMesh
from fealpy.mesh import HalfEdgeMesh2d, PolygonMesh
from fealpy.mesh import TriangleMesh
from mpl_toolkits.mplot3d import Axes3D

构造性方法

三角形网格

# 创建一个三角形网格
node = np.array([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)], dtype=np.float64)  # 结点， 容易知道我们的区域是[0,1]^2
cell = np.array([(1, 2, 0), (3, 0, 2)], dtype=np.int_)  # 单元, 我们可以画个图看一下代表的意思
mesh = TriangleMesh(node, cell)  #创建三角形网格

我们已经创建了一个三角形网格，为了更好地看到节点，网格的实际含义， 我们便画出图观察一下。 首先拿出网格的所有实体（网格，单元，边）。

node = mesh.entity('node')
edge = mesh.entity('edge')
cell = mesh.entity('cell')

然后利用 mesh 类中的函数 add_plot 画出相应的图

fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

mesh.find_node(axes, showindex=True)
mesh.find_edge(axes, showindex=True)
mesh.find_cell(axes, showindex=True)

从图中我们可以看到黄色的点为 cell 的编号， 0号三角形的编号是1->2->0, 1号三角形的编号是 0->3->2。
注意：cell 局部编号是逆时针排列， 对于边界边排列，我们要求左手边是边界的内部。
Fealpy 中的几种接口为

Fealpy 中构造了一个网格加密的函数 mesh.uniform_refine(n) ， 其核心想法便是把边的中点连接起来, n 指的是加密的次数， 比如：

mesh.uniform_refine(n=5)
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

我们可以在此基础上得到如下网格：

正方形网格

类似地， 我们可以利用QuadrangleMesh 来创建正方形网格：

node = np.array([(0, 0), (1, 0), (1, 1), (0, 1)], dtype=np.float64)
cell = np.array([(0, 1, 2, 3)], dtype=np.int_)

mesh =QuadrangleMesh(node, cell)
mesh.uniform_refine(2)
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

结果如下图：

注意到上面的网格是正方形的，那么我们能否创建长方形网格呢？ 答案是肯定的，我们不在使用原始的构造加密方式， 我们采用 mesh 的一个类来完成这件事情 ：Factory 。

利用 Mesh Factory 生成网格

from fealpy.mesh import MeshFactory
from matplotlib import pyplot as plt

mf = MeshFactory()

# 生成三角形网格
box = [0, 1, 0, 1]
mesh = mf.boxmesh2d(box, nx=4, ny=5, meshtype='tri')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

结果为

矩形网格

# 生成四边形网格
box = [0, 2, 0, 2]
mesh = mf.boxmesh2d(box, nx=4, ny=10, meshtype='quad')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

多项式网格

# 生成多项式网格

box = [0, 1, 0, 1]
mesh = mf.boxmesh2d(box, nx=4, ny=4, meshtype='poly')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

非结构三角网格

# 生成非结构的网格
box = [0, 1, 0, 1]
mesh = mf.triangle(box, h=0.1, meshtype='tri')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

非结构多项式网格

box = [0, 1, 0, 1]
mesh = mf.triangle(box, h=0.1, meshtype='poly')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

鱼骨型网格

# 鱼骨型
box = [0, 1, 0, 1]
mesh = mf.special_boxmesh2d(box, n=10, meshtype='fishbone')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

米字型网格

# 米字型
box = [0, 1, 0, 1]
mesh = mf.special_boxmesh2d(box, n=10, meshtype='rice')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

交叉型网格

# 交叉型
box = [0, 1, 0, 1]
mesh = mf.special_boxmesh2d(box, n=10, meshtype='cross')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

非一致网格

# 非一致型
box = [0, 1, 0, 1]
mesh = mf.special_boxmesh2d(box, n=10, meshtype='nonuniform')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

圆上的网格

# 圆上的网格
mesh = mf.unitcirclemesh(h=0.1, meshtype='tri')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

mesh = mf.unitcirclemesh(h=0.1, meshtype='poly')
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

L-shape 网格

# L-shape
mesh = mf.lshape_mesh(n=2)
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

Polygon 网格

mesh = mf.polygon_mesh()
mesh.uniform_refine(n=3)  # 网格加密次数
NC = mesh.number_of_cells()
print('Number of cells:', NC)
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)

四面体网格

# 四面体网格
box = [0, 1, 0, 1, 0, 1]
mesh = mf.boxmesh3d(box, nx=4, ny=4, nz=4, meshtype='tet')
fig = plt.figure()
axes = fig.gca(projection='3d')
mesh.add_plot(axes)

六面体网格

# 六面体网格
box = [0, 1, 0, 1, 0, 1]
mesh = mf.boxmesh3d(box, nx=4, ny=4, nz=4, meshtype='hex')
fig = plt.figure()
axes = fig.gca(projection='3d')
mesh.add_plot(axes)

利用Threshold构造你想要的网格

我们之前给出过L-shape的网格，但是有些时候想要变换 L-shape的方向，这时候便要借助 Threshold 函数来完成。

    # mesh = mf.boxmesh2d([-1, 1, -1, 1], nx=10, ny=10, meshtype='tri',
    #                    threshold=lambda p: (p[..., 0] > 0.0) & (p[..., 1] < 0.0))
    mesh = mf.boxmesh2d([-1, 1, -1, 1], nx=10, ny=10, meshtype='tri', threshold=lambda p: (p[..., 0]<0.0) & (p[..., 1]<0.0))
    mesh.add_plot(plt)

    mesh = mf.boxmesh2d([-1, 1, -1, 1], nx=100, ny=100, meshtype='tri', threshold=lambda p: p[..., 0] ** 2 + p[..., 1] ** 2 <0.5)
    mesh.add_plot(plt)

上图我们本意是想挖掉一个圆之后再离散，但这种写法显然是先离散，再将 threshold 中的区域去掉。 先挖掉一个圆再离散的话显然是自适应网格，怎么做我还不清楚。但是在内部挖掉一个正方形还是没有什么问题的。

总结

由上述例子我们可以看出， Fealpy 还是很强大的， 并且还在不断完善中。希望有更多的人使用Fealpy。下一节我们将介绍拉格朗日函数有限元空间的构造。