fealpy网格生成

构造性方法

import sys
import numpy as np
import matplotlib.pyplot as plt
from fealpy.mesh import TriangleMesh

node = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [1.0, 1.0],
    [0.0, 1.0]], dtype=np.float) # (NN, 2)

cell = np.array([[1, 2, 0], [3, 0, 2]], dtype=np.int) # (NC, 3)


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

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)
plt.show()

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

mesh.uniform_refine(n=1)
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)
plt.show()

正方形网格

import sys
import numpy as np
import matplotlib.pyplot as plt
from fealpy.mesh import QuadrangleMesh

node = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [1.0, 1.0],
    [0.0, 1.0]], dtype=np.float) 

cell = np.array([[0,1, 2,3]], dtype=np.int)


mesh =QuadrangleMesh(node, cell)

利用 Mesh Factory 生成网格

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

box = [0, 1, 0, 1]
mesh = MF.boxmesh2d(box, nx=2, ny=5, meshtype='tri')  #quad、poly

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)
plt.show()

非结构三角网格

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

box = [0, 1, 0, 1]
mesh = MF.triangle(box, h=0.2, meshtype='tri')

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)
plt.show()

在这里调用了MeshPy (tician.de), 其中安装需要c++ 14.0+函数库，[解决方案](pip错误“Microsoft Visual C++ 14.0 is required.”解决办法 - 知乎 (zhihu.com))

鱼骨型网格

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

box = [0, 1, 0, 1]
mesh = MF.special_boxmesh2d(box, n=4, meshtype='fishbone')

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)
plt.show()

米字型网格

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

box = [0, 1, 0, 1]
mesh = MF.special_boxmesh2d(box, n=4, meshtype='rice')

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)
plt.show()

交叉型网格

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

box = [0, 1, 0, 1]
mesh = MF.special_boxmesh2d(box, n=3, meshtype='cross')

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)
plt.show()

非一致网格

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

box = [0, 10, 0, 10]
mesh = MF.special_boxmesh2d(box, n=3, meshtype='nonuniform')

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)
plt.show()

圆上的网格

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

box = [0, 10, 0, 10]
mesh = MF.unitcirclemesh(h=0.4, meshtype='tri')
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)
plt.show()

L-shape 网格

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

box = [0, 10, 0, 10]
mesh = MF.lshape_mesh(n=1)
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)
plt.show()

四面体六面体网格

from fealpy.mesh import MeshFactory as MF
from matplotlib import pyplot as plt
box = [0, 1, 0, 1,0,1]
mesh = MF.boxmesh3d(box, nx=2, ny=2, nz=2, meshtype='tet') #hex
fig = plt.figure()
axes = fig.gca(projection='3d')
mesh.add_plot(axes)
plt.show()

挖掉一个圆网格

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

mesh = MF.boxmesh2d([-1, 1, -1, 1], nx=20, ny=20, meshtype='tri', threshold=lambda p: p[..., 0] ** 2 + p[..., 1] ** 2 <0.5)
mesh.add_plot(plt)
fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)