前言

安装fipy自带的库不能解决三维问题的显示和任意边界区域的设定。

• 前者需要额外安装mayavi
• 后者需要额外安装gmsh.

fipy

FiPy是基于标准有限体积（FV）方法的用Python编写的面向对象的偏微分方程（PDE）求解器。

安装

pip install future
pip install fipy

示例代码

from fipy import CellVariable, Grid2D, Viewer, TransientTerm, DiffusionTerm
from fipy.tools import numerix
from fipy import input

nx = 20
ny = nx
dx = 1.
dy = dx
L = dx * nx
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

phi = CellVariable(name = "solution variable",
                   mesh = mesh,
                   value = 0.)

D = 1.
eq = TransientTerm() == DiffusionTerm(coeff=D)                  

valueTopLeft = 0
valueBottomRight = 1

X, Y = mesh.faceCenters
facesTopLeft = ((mesh.facesLeft & (Y > L / 2))
                | (mesh.facesTop & (X < L / 2)))
facesBottomRight = ((mesh.facesRight & (Y < L / 2))
                    | (mesh.facesBottom & (X > L / 2)))



phi.constrain(valueTopLeft, facesTopLeft)
phi.constrain(valueBottomRight, facesBottomRight)

if __name__ == '__main__':
    viewer = Viewer(vars=phi, datamin=0., datamax=1.)
    viewer.plot()


timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)
steps = 10
from builtins import range
for step in range(steps):
    eq.solve(var=phi,
             dt=timeStepDuration)
    if __name__ == '__main__':
        viewer.plot()
    
if __name__ == '__main__':
    input("Implicit transient diffusion. Press <return> to proceed...")

mayavi

mayavi是三维问题的可视化显示界面的库。

安装

依赖VTK, traits,PyQt4,numpy.首先下载及安装PyQt4
需要从Python Extension Packages下载需要的版本安装.例如：

pip install PyQt4?4.11.4?cp36?cp36m?win32.whl

需要注意的是，PyQt4在笔者写作时没有python3.8及以上的版本。
之后输入以下代码安装mayavi及其依赖。

pip install mayavi

三维问题方体区域示例代码

from fipy import Grid3D, CellVariable, Viewer, DiffusionTerm, MayaviClient, TSVViewer
from fipy.tools import numerix
from fipy import input

# 网格参数
nx = 20
ny = 20
nz = 20 
dx = 1.
dy = dx
dz = dy 
L = nx * dx
# 三维网格
mesh = Grid3D(dx=dx, dy=dy, dz=dz, nx=nx,ny=ny,nz=nz)
# 单元格
phi = CellVariable(name = "solution variable",
                   mesh = mesh,
                   value = 0.)

# 所有方位 Back, Bottom, Down, Front, Left, Right, Top, Up

# 边值条件的值
valueFront  = 0
valueBack   = 0
valueLeft   = 0
valueRight  = 0
valueUp     = 1
valueDown   = 0

# 边界
X, Y, Z = mesh.faceCenters

facesFront  = ((mesh.facesFront ))
facesBack   = ((mesh.facesBack  ))
facesLeft   = ((mesh.facesLeft  ))
facesRight  = ((mesh.facesRight ))
facesUp     = ((mesh.facesUp    & ( Z < L/2)))
facesDown   = ((mesh.facesDown  ))

# 构造边值条件
phi.faceGrad.constrain(valueFront, facesFront)
phi.faceGrad.constrain(valueBack , facesBack )
phi.faceGrad.constrain(valueLeft , facesLeft )
phi.faceGrad.constrain(valueRight, facesRight)      # 第二边值
phi.constrain(valueUp   , facesUp   )               # 第一边值
phi.constrain(valueDown , facesDown )

# 观察器
viewer = MayaviClient(vars=phi,
                    limits={'ymin': 0.1, 'ymax': 0.9},
                    datamin=0, datamax=1,
                    title="MayaviClient test")
# 求稳态解
DiffusionTerm().solve(var=phi)      
# 保存数据
TSVViewer(vars=(phi, phi.grad)).plot(filename="myTSV.tsv")

# 绘出图像
viewer.plot()
viewer._promptForOpinion()

Gmsh

Gmsh是一个自动的三维有限元网格生成带有内置在CAD和后期处理器。

安装

安装gmsh

pip install gmsh

并从gmsh下载gmsh的压缩包。将解压后的文件的gmsh.exe放到python的Script文件夹下1。

二维示例代码（圆域）

cellSize = 0.05
radius = 1.

from fipy import CellVariable, Gmsh2D, TransientTerm, DiffusionTerm, Viewer
from fipy.tools import numerix

mesh = Gmsh2D('''
              cellSize = %(cellSize)g;
              radius = %(radius)g;
              Point(1) = {0, 0, 0, cellSize};
              Point(2) = {-radius, 0, 0, cellSize};
              Point(3) = {0, radius, 0, cellSize};
              Point(4) = {radius, 0, 0, cellSize};
              Point(5) = {0, -radius, 0, cellSize};
              Circle(6) = {2, 1, 3};
              Circle(7) = {3, 1, 4};
              Circle(8) = {4, 1, 5};
              Circle(9) = {5, 1, 2};
              Line Loop(10) = {6, 7, 8, 9};
              Plane Surface(11) = {10};
              ''' % locals()) 

phi = CellVariable(name = "solution variable",
                   mesh = mesh,
                   value = 0.) 

phi = CellVariable(name = "solution variable",
                   mesh = mesh,
                   value = 0.) 

viewer = None
from fipy import input
if __name__ == '__main__':
    try:
        viewer = Viewer(vars=phi, datamin=-1, datamax=1.)
        viewer.plotMesh()
        input("Irregular circular mesh. Press <return> to proceed...") 
    except:
        print("Unable to create a viewer for an irregular mesh (try Matplotlib2DViewer or MayaviViewer)")

D = 1.
eq = TransientTerm() == DiffusionTerm(coeff=D)

X, Y = mesh.faceCenters 
phi.constrain(X, mesh.exteriorFaces) 

timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
steps = 10
from builtins import range
for step in range(steps):
    eq.solve(var=phi,
             dt=timeStepDuration) 
    if viewer is not None:
        viewer.plot() 

三维示例代码（方域）

方域2

from fipy import CellVariable, Gmsh3D, TransientTerm, DiffusionTerm, Viewer, MayaviClient
from fipy.tools import numerix
from fipy import input

cellSize = 0.05
mesh = Gmsh3D('''
              cellSize = %(cellSize)g;
              Point(1) = {0, 0, 0, cellSize};
              Point(2) = {1, 0, 0, cellSize};
              Point(3) = {1, 1, 0, cellSize};
              Point(4) = {0, 1, 0, cellSize};
              Point(5) = {0, 0, 1, cellSize};
              Point(6) = {1, 0, 1, cellSize};
              Point(7) = {1, 1, 1, cellSize};
              Point(8) = {0, 1, 1, cellSize};
              Line(9)  = {1, 2};
              Line(10) = {2, 3};
              Line(11) = {3, 4};
              Line(12) = {4, 1};
              Line(13) = {5, 6};
              Line(14) = {6, 7};
              Line(15) = {7, 8};
              Line(16) = {8, 5};
              Line(17) = {1, 5};
              Line(18) = {2, 6};
              Line(19) = {3, 7};
              Line(20) = {4, 8};
              Line Loop(21) = {9, 10, 11, 12};
              Line Loop(22) = {13, 14, 15, 16};
              Line Loop(23) = {9, 18, -13, -17};
              Line Loop(24) = {10, 19, -14, -18};
              Line Loop(25) = {11, 20, -15, -19};
              Line Loop(26) = {12, 17, -16, -20};
              Plane Surface(27) = {21};
              Plane Surface(28) = {22};
              Plane Surface(29) = {23};
              Plane Surface(30) = {24};
              Plane Surface(31) = {25};
              Plane Surface(32) = {26};
              Surface Loop(33) = {27, 28, 29, 30, 31, 32};
              Volume(34) = {33};
              ''' % locals()) 

phi = CellVariable(name = "solution variable",
                   mesh = mesh,
                   value = 0.) 

viewer = None
if __name__ == '__main__':
    try:
        viewer = MayaviClient(vars=phi,
                    limits={'ymin': 0.1, 'ymax': 0.9},
                    datamin=0, datamax=1,
                    title="MayaviClient test")
        viewer.plotMesh()
        X, Y, Z = mesh.faceCenters 
        phi.constrain(X, mesh.exteriorFaces) 
        DiffusionTerm().solve(phi)
        viewer.plot()
        input("Irregular circular mesh. Press <return> to proceed...") 
    except:
        print("Unable to create a viewer for an irregular mesh (try Matplotlib2DViewer or MayaviViewer)")

感谢阅读！！！

多说一句，很多人学Python过程中会遇到各种烦恼问题，没有人解答容易放弃。小编是一名python开发工程师，这里有我自己整理了一套最新的python系统学习教程，包括从基础的python脚本到web开发、爬虫、数据分析、数据可视化、机器学习等。想要这些资料的可以关注小编，并在后台私信小编：“01”即可领取。